Paqueterias

Definimos los paquetes que se utilizarán

## Repositorio CRAN
cran_packages <- c("bookdown", "knitr", "tidyverse", "plyr", "grid",
                   "gridExtra", "kableExtra", "xtable", "ggpubr")
## Repositorio Bioconductor
bioc_packages <- c("phyloseq", "dada2", "DECIPHER", "phangorn", "ggpubr",
                   "BiocManager","DESeq2", "microbiome", "philr")
## Repositorio GitHub
git_source <- c("twbattaglia/btools", "gmteunisse/fantaxtic",
                "MadsAlbertsen/ampvis2", "opisthokonta/tsnemicrobiota")
# fuente/nombre del paquete
git_packages <- c("btools", "fantaxtic", "ampvis2", "tsnemicrobiota")

Instalamos los paquetes, para ello es necesario que descomentes las siguientes líneas.

# # Intalar paquetes CRAN
# .inst <- cran_packages %in% installed.packages()
# if(any(!.inst)) {
#   install.packages(cran_packages[!.inst])
# }
# # Intalar paquetes BioConductor
# if (!requireNamespace("BiocManager", quietly = TRUE))
#     install.packages("BiocManager")
# .inst <- bioc_packages %in% installed.packages()
# if(any(!.inst)) {
#   BiocManager::install(bioc_packages[!.inst])
# }
# # Instalar paquetes GitHub
# .inst <- git_source %in% installed.packages()
# install.packages("devtools")
# if(any(!.inst)) {
#   devtools::install_github(git_source[!.inst])
# }
# library("devtools")
# install_github("opisthokonta/tsnemicrobiota")
# install.packages("remotes")
# remotes::install_github("kasperskytte/ampvis2")
# devtools::install_github('twbattaglia/btools')

Cargamos los paquetes

sapply(c(cran_packages, bioc_packages, git_packages), require, 
       character.only = TRUE)
##       bookdown          knitr      tidyverse           plyr           grid 
##           TRUE           TRUE           TRUE           TRUE           TRUE 
##      gridExtra     kableExtra         xtable         ggpubr       phyloseq 
##           TRUE           TRUE           TRUE           TRUE           TRUE 
##          dada2       DECIPHER       phangorn         ggpubr    BiocManager 
##           TRUE           TRUE           TRUE           TRUE           TRUE 
##         DESeq2     microbiome          philr         btools      fantaxtic 
##           TRUE           TRUE           TRUE           TRUE           TRUE 
##        ampvis2 tsnemicrobiota 
##           TRUE           TRUE
library("phyloseq")
library("ggplot2")
library("readr")
library("patchwork")
library("pheatmap")
library("microbiome")
library("xtable")
library("tsnemicrobiota")
library("Rtsne")
library("ampvis2")
library("btools")
library("RColorBrewer")
library("igraph")
library("vegan")
library("factoextra")
library("pbkrtest")
library("BiodiversityR")

Directorio de trabajo

Ajustamos directorio de trabajo donde se encuentran las lecturas “reads”

# IMPORTANTE
# Ajustamos path de trabajo según tu PC
# Mostramos directorio actual
getwd()
## [1] "C:/Users/alber/OneDrive - CINVESTAV/Documentos/Doctorado/ProyectosR/Calakmul/Shotgun"
# Si es necesario, cambiar la ruta donde están los archivos almacenados. 
# "." significa el directorio actual. 
workingDir <-  "."
setwd(workingDir)
# Creamos directorio de Resultados y otras capertas
dir.create("Resultados/")
## Warning in dir.create("Resultados/"): 'Resultados' already exists
dir.create("Resultados/Abundancias/")
## Warning in dir.create("Resultados/Abundancias/"): 'Resultados\Abundancias'
## already exists
dir.create("Resultados/AbundanciasAV2/")
## Warning in dir.create("Resultados/AbundanciasAV2/"): 'Resultados\AbundanciasAV2'
## already exists
dir.create("Resultados/BetaDiversidad/")
## Warning in dir.create("Resultados/BetaDiversidad/"): 'Resultados\BetaDiversidad'
## already exists
# Asignamos la carpeta donde están las lecturas crudas
miseq_path <- "RawData"
list.files(miseq_path)
## character(0)

Importación y manipulación de archivo kraken-biom

Metemos a la variable merged_metagenomes el resultado de kraken-biom, esto resulta en un objeto phyloseq

merged_metagenomes <- import_biom("CalakmulShotGun.biom")
class(merged_metagenomes)
## [1] "phyloseq"
## attr(,"package")
## [1] "phyloseq"
head(merged_metagenomes@tax_table@.Data)
##         Rank1         Rank2               Rank3                   
## 356     "k__Bacteria" "p__Proteobacteria" "c__Alphaproteobacteria"
## 41294   "k__Bacteria" "p__Proteobacteria" "c__Alphaproteobacteria"
## 374     "k__Bacteria" "p__Proteobacteria" "c__Alphaproteobacteria"
## 1437360 "k__Bacteria" "p__Proteobacteria" "c__Alphaproteobacteria"
## 1274631 "k__Bacteria" "p__Proteobacteria" "c__Alphaproteobacteria"
## 2057741 "k__Bacteria" "p__Proteobacteria" "c__Alphaproteobacteria"
##         Rank4            Rank5                  Rank6              
## 356     "o__Rhizobiales" "f__"                  "g__"              
## 41294   "o__Rhizobiales" "f__Bradyrhizobiaceae" "g__"              
## 374     "o__Rhizobiales" "f__Bradyrhizobiaceae" "g__Bradyrhizobium"
## 1437360 "o__Rhizobiales" "f__Bradyrhizobiaceae" "g__Bradyrhizobium"
## 1274631 "o__Rhizobiales" "f__Bradyrhizobiaceae" "g__Bradyrhizobium"
## 2057741 "o__Rhizobiales" "f__Bradyrhizobiaceae" "g__Bradyrhizobium"
##         Rank7            
## 356     "s__"            
## 41294   "s__"            
## 374     "s__"            
## 1437360 "s__erythrophlei"
## 1274631 "s__icense"      
## 2057741 "s__sp. SK17"

En los datos, el nombre como tal empieza en la letra numero 4, tiene 3 antes, por ejemplo “k__“, con el siguiente comando eliminamos estas 3 primeras letras y empieza en la 4ta.

merged_metagenomes@tax_table@.Data <- substring(merged_metagenomes@tax_table@.Data, 4)
# Cambiamos el nombre de las columnas con el siguiente comando
colnames(merged_metagenomes@tax_table@.Data) <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species")
head(merged_metagenomes@tax_table@.Data)
##         Kingdom    Phylum           Class                 Order        
## 356     "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
## 41294   "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
## 374     "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
## 1437360 "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
## 1274631 "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
## 2057741 "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
##         Family              Genus            Species       
## 356     ""                  ""               ""            
## 41294   "Bradyrhizobiaceae" ""               ""            
## 374     "Bradyrhizobiaceae" "Bradyrhizobium" ""            
## 1437360 "Bradyrhizobiaceae" "Bradyrhizobium" "erythrophlei"
## 1274631 "Bradyrhizobiaceae" "Bradyrhizobium" "icense"      
## 2057741 "Bradyrhizobiaceae" "Bradyrhizobium" "sp. SK17"
# Lo metemos en una tabla
TablaMergedMetagenomes <- merged_metagenomes@tax_table@.Data
# Guardamos
write.table(TablaMergedMetagenomes, file = "CalakmulTablaMergedMetagenomes.txt", 
            row.names=TRUE, col.names=NA, quote=FALSE, sep="\t")

Podemos ver los únicos Kingdom, Phylum, Class, etc…

head(unique(merged_metagenomes@tax_table@.Data[,"Kingdom"]))
## [1] "Bacteria"
head(unique(merged_metagenomes@tax_table@.Data[,"Phylum"]))
## [1] "Proteobacteria"      "Actinobacteria"      "Firmicutes"         
## [4] "Deinococcus-Thermus" "Cyanobacteria"       "Chloroflexi"
head(unique(merged_metagenomes@tax_table@.Data[,"Class"]))
## [1] "Alphaproteobacteria"   "Betaproteobacteria"    "Gammaproteobacteria"  
## [4] "Deltaproteobacteria"   "Epsilonproteobacteria" "Acidithiobacillia"
head(unique(merged_metagenomes@tax_table@.Data[,"Order"]))
## [1] "Rhizobiales"      "Sphingomonadales" "Rhodospirillales" "Rhodobacterales" 
## [5] "Caulobacterales"  "Rickettsiales"
head(unique(merged_metagenomes@tax_table@.Data[,"Family"]))
## [1] ""                    "Bradyrhizobiaceae"   "Rhizobiaceae"       
## [4] "Methylobacteriaceae" "Hyphomicrobiaceae"   "Phyllobacteriaceae"
head(unique(merged_metagenomes@tax_table@.Data[,"Genus"]))
## [1] ""                 "Bradyrhizobium"   "Rhodopseudomonas" "Bosea"           
## [5] "Afipia"           "Nitrobacter"
head(unique(merged_metagenomes@tax_table@.Data[,"Species"]))
## [1] ""               "erythrophlei"   "icense"         "sp. SK17"      
## [5] "lablabi"        "diazoefficiens"

Como en Kingdom hay “solo” bacterias, nos vamos a quedar solo con eso, y así quitar lo que no entra en esta categoria.

merged_metagenomes <- subset_taxa(merged_metagenomes, Kingdom == "Bacteria")

Para ver cuantos de un cierto taxón hay dentro de otro

# Aquí vemos cuantos Firmicutes hay en el Phylum
sum(merged_metagenomes@tax_table@.Data[,"Phylum"] == "Firmicutes")
## [1] 708

Añadimos metadatos

MetaShotgun <- data.frame(SAMPLE=c("Ag1", "Ag2", "Ag3"), 
                          SITE=c("2", "3", "1"), 
                          pH=c(6.2, 6.7, 6.1), 
                          N_TOTAL=c(0.36, 0.22, 0.29),
                          P_OLSEN=c(29.6, 13.5, 4.0),
                          ORGANIC=c(12.0, 13.5, 24.2),
                          SITE_ID=c("Ag1", "Ag2", "Ag3"))
rownames(MetaShotgun) <- MetaShotgun$SAMPLE
rownames(merged_metagenomes) %in% rownames(MetaShotgun)
## logical(0)
all(rownames(merged_metagenomes) %in% MetaShotgun$SAMPLE)
## [1] TRUE
merged_metagenomes@sam_data <- sample_data(MetaShotgun)

Ensamble y manipulación de objetos de Phyloseq

Resumen del objeto Phyloseq

# Nos dice el número de taxas en las muestras y en cuantos rangos taxonómicos
merged_metagenomes 
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 4457 taxa and 3 samples ]
## sample_data() Sample Data:       [ 3 samples by 7 sample variables ]
## tax_table()   Taxonomy Table:    [ 4457 taxa by 7 taxonomic ranks ]
# Cantidad de elementos en cada muestra
sample_sums(merged_metagenomes)
##     Ag1     Ag2     Ag3 
## 1633465 1844772 1611481
# Un resumen de cada muestra, indicando cuantos mínimos medios y máximos hay para cada taxa
summary(merged_metagenomes@otu_table@.Data) 
##       Ag1               Ag2               Ag3         
##  Min.   :    0.0   Min.   :    0.0   Min.   :    0.0  
##  1st Qu.:    8.0   1st Qu.:    9.0   1st Qu.:    8.0  
##  Median :   39.0   Median :   38.0   Median :   37.0  
##  Mean   :  366.5   Mean   :  413.9   Mean   :  361.6  
##  3rd Qu.:  307.0   3rd Qu.:  328.0   3rd Qu.:  295.0  
##  Max.   :40397.0   Max.   :49917.0   Max.   :39700.0

Gráfico

# Generamos un data.frame con datos del merge_metagenomes
deep <- data.frame(Samples = sample_names(merged_metagenomes),
                    Reads = sample_sums(merged_metagenomes)/1000000)
deep
##     Samples    Reads
## Ag1     Ag1 1.633465
## Ag2     Ag2 1.844772
## Ag3     Ag3 1.611481
# Graficamos
COLOR = c("#a6cee3", "#b2df8a", "#fdbf6f")
plot.reads <- ggplot(data = deep, mapping = aes(x = Samples, y = Reads, fill = Samples)) +
  theme_classic() +
  labs(x="Samples", y="Million reads",
       title = "Total readings of analyzed samples") +
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold")) +
  geom_col(position = position_dodge(0.9), width = 0.9) +
  scale_fill_brewer(palette = "Set1")
# Verla en el markdown
plot.reads

pdf("PlotReads.pdf", width=10, height=7)
plot.reads
dev.off()
## png 
##   2
svg("PlotReads.svg", width=10, height=7)
plot.reads
dev.off()
## png 
##   2

Estructura y manipulación de un objeto phyloseq

Muchas veces queremos analizar un sub conjunto de las muestras en nuestro objeto phyloseq, o bien, queremos seleccionar ciertos grupos taxonómicos para análisis posteriores. Phyloseq nos permite hacer todo tipo de filtros para llevar esto a cabo. Veamos dónde se almacena la información en phyloseq.

# Número de taxas
ntaxa(merged_metagenomes)
## [1] 4457
# Número de muestras
nsamples(merged_metagenomes)
## [1] 3
# Nombre de las muestras
sample_names(merged_metagenomes)
## [1] "Ag1" "Ag2" "Ag3"
# Niveles taxonómicos 
rank_names(merged_metagenomes)
## [1] "Kingdom" "Phylum"  "Class"   "Order"   "Family"  "Genus"   "Species"
# Nombre de las variables de los metadatos
sample_variables(merged_metagenomes)
## [1] "SAMPLE"  "SITE"    "pH"      "N_TOTAL" "P_OLSEN" "ORGANIC" "SITE_ID"

También podemos ver un resumen estadístico de la tabla de OTUs

# Resumen estadístico
summary(merged_metagenomes@otu_table@.Data)
##       Ag1               Ag2               Ag3         
##  Min.   :    0.0   Min.   :    0.0   Min.   :    0.0  
##  1st Qu.:    8.0   1st Qu.:    9.0   1st Qu.:    8.0  
##  Median :   39.0   Median :   38.0   Median :   37.0  
##  Mean   :  366.5   Mean   :  413.9   Mean   :  361.6  
##  3rd Qu.:  307.0   3rd Qu.:  328.0   3rd Qu.:  295.0  
##  Max.   :40397.0   Max.   :49917.0   Max.   :39700.0

La tabla de taxonomía relaciona el nombre de las taxa con el linaje taxonómico de éstas, por ejemplo, vincula una variante de secuencia, o ASV, con los rangos taxonómicos desde Reino hasta Género o Especie dependiendo del nivel de resolución del análisis.

# Tabla de taxonomía
tax_table(merged_metagenomes)[1:5, 1:4]
## Taxonomy Table:     [5 taxa by 4 taxonomic ranks]:
##         Kingdom    Phylum           Class                 Order        
## 356     "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
## 41294   "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
## 374     "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
## 1437360 "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"
## 1274631 "Bacteria" "Proteobacteria" "Alphaproteobacteria" "Rhizobiales"

Ahora, el objeto phyloseq se ha vuelto un estándar en la industria ya que otros paquetes ahora usan esta estructura de datos para sus propias funciones. Uno de esos paquetes es microbiome y ampvis. Podemos fácilmente obtener un resumen global de nuestro objeto phyloseq usando la función summarize_phyloseq.

summarize_phyloseq(merged_metagenomes)
## [[1]]
## [1] "1] Min. number of reads = 1611481"
## 
## [[2]]
## [1] "2] Max. number of reads = 1844772"
## 
## [[3]]
## [1] "3] Total number of reads = 5089718"
## 
## [[4]]
## [1] "4] Average number of reads = 1696572.66666667"
## 
## [[5]]
## [1] "5] Median number of reads = 1633465"
## 
## [[6]]
## [1] "7] Sparsity = 0.0270735173135891"
## 
## [[7]]
## [1] "6] Any OTU sum to 1 or less? YES"
## 
## [[8]]
## [1] "8] Number of singletons = 65"
## 
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n        (i.e. exactly one read detected across all samples)1.4583800762845"
## 
## [[10]]
## [1] "10] Number of sample variables are: 7"
## 
## [[11]]
## [1] "SAMPLE"  "SITE"    "pH"      "N_TOTAL" "P_OLSEN" "ORGANIC" "SITE_ID"

Este comando nos muestra el mínimo y máximo de reads, número total y promedio de reads, etc. También muestra los encabezados de las columans en la tabla de metadata. Veamos ahora una tabla que mezcla metadata, taxonomía y abundancia del taxon más abundante de cada muestra.

df <- psmelt(merged_metagenomes)
head(df)
##          OTU Sample Abundance SAMPLE SITE  pH N_TOTAL P_OLSEN ORGANIC SITE_ID
## 4800    1883    Ag2     49917    Ag2    3 6.7    0.22    13.5    13.5     Ag2
## 4798    1883    Ag1     40397    Ag1    2 6.2    0.36    29.6    12.0     Ag1
## 4799    1883    Ag3     39700    Ag3    1 6.1    0.29     4.0    24.2     Ag3
## 5126  191495    Ag2     36610    Ag2    3 6.7    0.22    13.5    13.5     Ag2
## 9155     374    Ag2     25994    Ag2    3 6.7    0.22    13.5    13.5     Ag2
## 11663 674703    Ag2     25729    Ag2    3 6.7    0.22    13.5    13.5     Ag2
##        Kingdom         Phylum               Class               Order
## 4800  Bacteria Actinobacteria      Actinobacteria    Streptomycetales
## 4798  Bacteria Actinobacteria      Actinobacteria    Streptomycetales
## 4799  Bacteria Actinobacteria      Actinobacteria    Streptomycetales
## 5126  Bacteria Actinobacteria     Thermoleophilia Solirubrobacterales
## 9155  Bacteria Proteobacteria Alphaproteobacteria         Rhizobiales
## 11663 Bacteria Proteobacteria Alphaproteobacteria         Rhizobiales
##                  Family          Genus       Species
## 4800  Streptomycetaceae   Streptomyces              
## 4798  Streptomycetaceae   Streptomyces              
## 4799  Streptomycetaceae   Streptomyces              
## 5126  Conexibacteraceae   Conexibacter        woesei
## 9155  Bradyrhizobiaceae Bradyrhizobium              
## 11663 Hyphomicrobiaceae    Rhodoplanes sp. Z2-YC6860

Ahora veamos cómo podemos filtrar y hacer subsetting de un objeto phyloseq. Esto lo hacemos con tres grupos de funciones, por ejemplo, filter, subset, y prune. + Filtrar se refiere a filtrar según alguna regla lógica. Ya lo hicimos en la parte de control de calidad cuando llamamos la función filter_taxa(psd1, function(x) mean(x) > 1e-5, TRUE). Acá le pedíamos a la función filter_taxa sobre el objeto merged_metagenomes, calculara la media de cuentas de lecturas para cada taxa y si este resultado era menor que 1e-5, lo eliminara. Ahora filtramos según abundancia. Primero transformamos en abundancia relativa y luego filtramos.

# Transformamos las cuentas en porcentaje
Porcentaje <- transform_sample_counts(merged_metagenomes, function(x) x / sum(x) * 100 )
# Filtramos las taxa con una abundancia inferior al 1%
Porcentaje.filtrado <- filter_taxa(Porcentaje, function(x) sum(x) > 1, TRUE)

¿Cuántas taxa permanecen en nuestro objeto phyloseq? Con una operación tan simple como la que acabamos de aplicar, nos damos cuenta que la mayoría de las taxa presentes en nuestras están en muy baja abundancia. Ahora imaginemos la situación donde queremos filtrar nuestro objeto pero en función de un taxon en específico.

# Ahora filtramos de acuerdo a Acidobacteriota
Porcentaje.filtrado.Proteobacteria <- subset_taxa(Porcentaje.filtrado, Phylum == "Proteobacteria")
# También podríamos todo lo que NO es Acidobacteriota
Porcentaje.filtrado.NoAcidobacteriota <- subset_taxa(Porcentaje.filtrado, Phylum != "Acidobacteriota")

Otra manera de filtrar un objeto phyloseq es en base a algún atributo presente en sample_data. Por ejemplo, con estos datos uno podría querer estudiar el microbioma de área por separado. Para esto crearíamos tres objetos phyloseq a partir de merged_metagenomes.

merged_metagenomes.Ag1 <- subset_samples(merged_metagenomes, SITE_ID == "Ag1")
merged_metagenomes.Ag2 <- subset_samples(merged_metagenomes, SITE_ID == "Ag2")
merged_metagenomes.Ag3 <- subset_samples(merged_metagenomes, SITE_ID == "Ag3")

Alternativamente, podríamos decidir estudiar solo dos de las 3 áreas

merged_metagenomes.Ag1_Ag2 <- subset_samples(subset_samples(merged_metagenomes, SITE_ID != "Ag3"))

Transformación y manipulación de datos

Siguimos viendo el objeto phyloseq

# Con esto nos dice cuántos datos vacios hay, en este caso serán todos aquellos que aparezcan como TRUE
summary(merged_metagenomes@tax_table@.Data== "") # Hay 1337 NAs
##   Kingdom          Phylum          Class           Order        
##  Mode :logical   Mode :logical   Mode :logical   Mode :logical  
##  FALSE:4457      FALSE:4457      FALSE:4313      FALSE:4445     
##                                  TRUE :144       TRUE :12       
##    Family          Genus          Species       
##  Mode :logical   Mode :logical   Mode :logical  
##  FALSE:4381      FALSE:4217      FALSE:3807     
##  TRUE :76        TRUE :240       TRUE :650
#Actualizamos la variable merged_metagenomes indicando que se incluya todo lo diferente (!=) a vacio "", para la variable "Genus"
merged_metagenomes <- subset_taxa(merged_metagenomes, Genus != "")
summary(merged_metagenomes@tax_table@.Data== "") # Ya no hay ninguno
##   Kingdom          Phylum          Class           Order        
##  Mode :logical   Mode :logical   Mode :logical   Mode :logical  
##  FALSE:4217      FALSE:4217      FALSE:4088      FALSE:4210     
##                                  TRUE :129       TRUE :7        
##    Family          Genus          Species       
##  Mode :logical   Mode :logical   Mode :logical  
##  FALSE:4206      FALSE:4217      FALSE:3791     
##  TRUE :11                        TRUE :426

Datos Absolutos_S y Relativos_S

Generamos una tabla de datos Absolutos_S y Relativos_S a nivel de “Phylum”, para posteriores análisis.

Relativos_S <- transform_sample_counts(merged_metagenomes, function(x) x*100 / sum(x))
# Valores Relativos_S
Glom_Rel_S <- tax_glom(physeq = Relativos_S, taxrank = 'Phylum')
# psmelt permite generar un data.frame a partir de un objeto de phyloseq
Relativos_S <- psmelt(Glom_Rel_S)
# Valores Absolutos_S
Glom_Abs_S <- tax_glom(physeq = merged_metagenomes, taxrank = "Phylum")
Absolutos_S <- psmelt(Glom_Abs_S)

Gráficas de valores Absolutos_S y Relativos_S

En eje X usamos la variable “SITE_ID”, y en Y es “Abundance” y rellenamos por “Phylum”

Abs.plot <- ggplot(data=Absolutos_S, aes(x=SITE_ID, y=Abundance, fill=Phylum)) + 
  geom_bar(aes(), stat="identity", position="stack") + 
  theme_classic() +
  labs(x="Samples", y="Absolute abundance",
       title = "Absolute abundance of analyzed samples") +
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold"))

Rel.plot <- ggplot(data=Relativos_S, aes(x=SITE_ID, y=Abundance, fill=Phylum)) + 
  geom_bar(aes(), stat="identity", position="fill") + 
  theme_classic() +
  labs(x="Samples", y="Relative abundance",
       title = "Relative abundance of analyzed samples") +
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold"))

Abs.plot | Rel.plot
Abundancia absoluta y relativa a nivel de *Phylum*.

Abundancia absoluta y relativa a nivel de Phylum.

pdf("Resultados/Abundancias/PlotRaw.pdf", width=10, height=7)
Abs.plot
dev.off()
## png 
##   2
svg("Resultados/Abundancias/PlotRaw.svg", width=10, height=7)
Abs.plot
dev.off()
## png 
##   2
pdf("Resultados/Abundancias/PlotRel.pdf", width=10, height=7)
Rel.plot
dev.off()
## png 
##   2
svg("Resultados/Abundancias/PlotRel.svg", width=10, height=7)
Rel.plot
dev.off()
## png 
##   2
pdf("Resultados/Abundancias/PlotRaw-Rel.pdf", width=15, height=7)
Abs.plot | Rel.plot
dev.off()
## png 
##   2
svg("Resultados/Abundancias/PlotRaw-Rel.svg", width=15, height=7)
Abs.plot | Rel.plot
dev.off()
## png 
##   2

Podemos hacerlo para las diferentes variables en X, Y y fill

names(Relativos_S)
##  [1] "OTU"       "Sample"    "Abundance" "SAMPLE"    "SITE"      "pH"       
##  [7] "N_TOTAL"   "P_OLSEN"   "ORGANIC"   "SITE_ID"   "Kingdom"   "Phylum"
names(Absolutos_S)
##  [1] "OTU"       "Sample"    "Abundance" "SAMPLE"    "SITE"      "pH"       
##  [7] "N_TOTAL"   "P_OLSEN"   "ORGANIC"   "SITE_ID"   "Kingdom"   "Phylum"

Renombramos aquellos “Phylum” que presenten una abundancia menor a 1 %.

Relativos_S$Phylum <- as.character(Relativos_S$Phylum) # Primero lo convertimos a letra
Relativos_S$Phylum[Relativos_S$Abundance < 1.0] <- "Phylum < 1.0% abund." # Renombramos
unique(Relativos_S$Phylum)
## [1] "Proteobacteria"       "Actinobacteria"       "Firmicutes"          
## [4] "Planctomycetes"       "Acidobacteria"        "Phylum < 1.0% abund."

Volvemos a graficar usando ese filtro de “Phylum < 1.0% abund.”

Rel.plot2 <- ggplot(data=Relativos_S, aes(x=SITE_ID, y=Abundance, fill=Phylum)) + 
  geom_bar(aes(), stat="identity", position="fill") + 
  theme_classic() +
  labs(x="Samples", y="Relative abundance",
       title = "Relative abundance of analyzed samples") +
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold"))

Abs.plot | Rel.plot2
Abundancia absoluta y relativa a nivel de *Phylum*. con Phylum < 1.0 % abundancia.

Abundancia absoluta y relativa a nivel de Phylum. con Phylum < 1.0 % abundancia.

pdf("Resultados/Abundancias/PlotRel2.pdf", width=10, height=7)
Rel.plot2
dev.off()
## png 
##   2
svg("Resultados/Abundancias/PlotRel2.svg", width=10, height=7)
Rel.plot2
dev.off()
## png 
##   2
pdf("Resultados/Abundancias/PlotRaw-Rel2.pdf", width=15, height=7)
Abs.plot | Rel.plot2
dev.off()
## png 
##   2
svg("Resultados/Abundancias/PlotRaw-Rel2.svg", width=15, height=7)
Abs.plot | Rel.plot2
dev.off()
## png 
##   2

Abundancias

Graficamos con ggplot2 para ver el comportamiento de las lecturas en cada una de las muestras

#Generamos un data.frame con datos del merge_metagenomes
deep <- data.frame(Samples = merged_metagenomes@sam_data@.Data[[1]],
                   Reads = sample_sums(merged_metagenomes))
deep
##     Samples   Reads
## Ag1     Ag1 1523762
## Ag2     Ag2 1726131
## Ag3     Ag3 1503461
#Graficamos
plot.reads <- ggplot(data = deep, mapping = aes(x = Samples, y = Reads, fill = Samples)) +
  theme_classic() +
  labs(x="Samples", y="Reads",
       title = "Total readings of analyzed samples") +
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
        plot.title = element_text(hjust = 0.5, size=14, face="bold")) +
  geom_col(position = position_dodge(0.9), width = 0.9) +
  scale_fill_brewer(palette = "Paired")
plot.reads
Distribución del total de lecturas en cada una de las muestras.

Distribución del total de lecturas en cada una de las muestras.

pdf("Resultados/Abundancias/PlotReads.pdf", width=10, height=7)
plot.reads
dev.off()
## png 
##   2
svg("Resultados/Abundancias/PlotReads.svg", width=10, height=7)
plot.reads
dev.off()
## png 
##   2

Nos quedamos con el nombre de aquellos Phylum con procentaje mayor a 1

OnePercentage <- unique(Relativos_S$Phylum[Relativos_S$Abundance > 1.0])
OnePercentage
## [1] "Proteobacteria" "Actinobacteria" "Firmicutes"     "Planctomycetes"
## [5] "Acidobacteria"

Abundancias a nivel de Phylum nuevamente para los de abundancia mayor al 1 %

for (i in OnePercentage){
  tryCatch({
    # Prueba 
    # i <- "Actinobacteriota"
    phy <- subset_taxa(merged_metagenomes, Kingdom == "Bacteria") # Nos quedamos solo con Bacterias
    phy <- subset_taxa(merged_metagenomes, Phylum == i) # Nos quedamos con el phylum de interés
    # Comprobamos
    paste("El siguiente gráfico es del Phylum", unique(phy@tax_table@.Data[,2])) 
    # Lo hacemos valores Relativos_S
    phy <- transform_sample_counts(phy, function(x) x*100 / sum(x) ) 
    # Cortamos a nivel de género
    phy <- tax_glom(phy, taxrank = "Genus")
    # Lo volvemos data.frame
    phy.t <- psmelt(phy)
    # Renombramos los géneros que tienen abundancia menor a 1 %
    phy.t$Genus[phy.t$Abundance < 1.0] <- "Genus < 1.0 % abundance"
    }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")}) 
    # Graficamos resécto a "SITE_ID"
    plot <- ggplot(data = phy.t, aes(x=SITE_ID, y=Abundance, fill=Genus)) + 
      geom_bar(aes(), stat="identity", position="fill") +
      theme(legend.position = "right", 
            legend.text = element_text(size=12, face="bold"),
            legend.title = element_text(size = 12, face = "bold"),
            axis.text = element_text(color = "black", size = 12, face = "bold"),
            axis.title = element_text(color = "black", size = 12, face = "bold"),
            plot.title = element_text(hjust = 0.5, size=14, face="bold")) +
      labs(x="Samples", y="Relative abundance",
          title = paste("Abundancia relativa dentro del phylum", i))
    print(plot)
    pdf(paste0("Resultados/Abundancias/PlotRel", i, ".pdf"), width = 30, height = 10)
    print(plot)
    dev.off()
    svg(paste0("Resultados/Abundancias/PlotRel", i, ".svg"), width = 30, height = 10)
    print(plot)
    dev.off()
}
Abundancia relativa de todos los *Phylum* presentes en mayor a 1 % de Abundancia.

Abundancia relativa de todos los Phylum presentes en mayor a 1 % de Abundancia.

Abundancia relativa de todos los *Phylum* presentes en mayor a 1 % de Abundancia.

Abundancia relativa de todos los Phylum presentes en mayor a 1 % de Abundancia.

Abundancia relativa de todos los *Phylum* presentes en mayor a 1 % de Abundancia.

Abundancia relativa de todos los Phylum presentes en mayor a 1 % de Abundancia.

Abundancia relativa de todos los *Phylum* presentes en mayor a 1 % de Abundancia.

Abundancia relativa de todos los Phylum presentes en mayor a 1 % de Abundancia.

Abundancia relativa de todos los *Phylum* presentes en mayor a 1 % de Abundancia.

Abundancia relativa de todos los Phylum presentes en mayor a 1 % de Abundancia.

Análisis de diversidad

¿Qué entendemos por diversidad? Al menos podemos conceptualizar diversidad a dos niveles: diversidad genética o morfológica, y biodiversidad. En el estudio de comunidades, tomamos prestado el concepto de biodiversidad de ecología de comunidades donde estamos interesados en la riqueza de especies (número de especies diferentes en una comunidad o diversidad alfa), en las diferencias y similitudes entre comunidades (diversidad beta), y en algunos casos en la diversidad total de un región o paisaje ecológico (landscape; diversidad gama). Medidas de riqueza, uniformidad, dominancia, diversidad filogenética (diversidad alfa)

Alfa diversidad

En el contexto metagenómico, medimos diversidad alfa \(\alpha\) usando una serie de medidas prestadas de ecología que nos permiten caracterizar una comunidad microbiana. phyloseq tiene una función muy útil que nos permite calcular y graficar hasta siete medidas, i .e., Observed (simplemente el número de taxa o riqueza), Chao1 (la riqueza ajustada por probabilidad de no observar especies), ACE (riqueza que toma en cuenta la abundancia relativa), Shannon (abundancia relativa de taxa), Simpson (1 - la probabilidad de que observemos aleatoriamente dos bacterias en una comunidad y que pertenezcan a diferentes especies ), Inverse Simpson ( 1 / Simpson), y Fisher (riqueza tomando en cuenta abundancia). En phyloseq simplemente llamamos la función plot_richness y podemos visualizar las medidas de diversidad.

alfa.diver <- plot_richness(physeq = merged_metagenomes, 
              color =  "SITE_ID",
              x = "SITE_ID",
              measures = c("Observed", "Chao1", "ACE", "Shannon", "Simpson", "Fisher"),
              title = "Alpha diversity indices for the Calakmul samples obtained by ShotGun") + 
  labs(x="Samples", y="Alpha Diversity Measure") +
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold")) +
  geom_boxplot(aes(fill = SITE_ID), alpha=.7) + 
  scale_color_brewer(palette = "Set1") + 
  scale_fill_brewer(palette = "Set1")
alfa.diver
Índices de diversidad alpha $lpha$.

Índices de diversidad alpha \(lpha\).

pdf("Resultados/PlotAlfa.pdf", width=20, height=10)
alfa.diver
dev.off()
## png 
##   2
svg("Resultados/PlotAlfa.svg", width=20, height=10)
alfa.diver
dev.off()
## png 
##   2

El paquete microbiome ofrece otras herramientas para evaluar diversidad que son accesibles fácilmente a través de su función global.

Alpha <- alpha(merged_metagenomes, index = "all")
head(Alpha)
##     observed    chao1 diversity_inverse_simpson diversity_gini_simpson
## Ag1     4122 4214.560                  352.4709              0.9971629
## Ag2     4131 4214.379                  306.7775              0.9967403
## Ag3     4100 4189.446                  345.0244              0.9971017
##     diversity_shannon diversity_fisher diversity_coverage evenness_camargo
## Ag1          6.885202         515.8116                231        0.2035280
## Ag2          6.811040         508.0487                212        0.1921895
## Ag3          6.869876         513.6507                228        0.2018350
##     evenness_pielou evenness_simpson evenness_evar evenness_bulla dominance_dbp
## Ag1       0.8271414       0.08550968     0.1350859      0.3609202    0.02651136
## Ag2       0.8180177       0.07426228     0.1312996      0.3506728    0.02891843
## Ag3       0.8258311       0.08415229     0.1360249      0.3592546    0.02640574
##     dominance_dmn dominance_absolute dominance_relative dominance_simpson
## Ag1    0.04183593              40397         0.02651136       0.002837114
## Ag2    0.05012771              49917         0.02891843       0.003259692
## Ag3    0.04186274              39700         0.02640574       0.002898346
##     dominance_core_abundance dominance_gini rarity_log_modulo_skewness
## Ag1                0.2763476      0.8010570                   2.054063
## Ag2                0.2906512      0.8116390                   2.052640
## Ag3                0.2796986      0.8035928                   2.053428
##     rarity_low_abundance rarity_rare_abundance
## Ag1            0.7161781             0.7006737
## Ag2            0.7008865             0.6776919
## Ag3            0.7180931             0.6992699

La función alpha nos da 22 medidas de diversidad que nos ayudan a entender la estructura de las comunidades microbianas. En general, estas medidas se dividen en riqueza, diversidad, dominancia, rareza, cobertura y uniformidad. El paquete microbiome ofrece funciones para calcular cada uno de estos aspectos de las comunidades microbianas.

# Riqueza
Riqueza <- richness(merged_metagenomes)
Riqueza
##     observed    chao1
## Ag1     4122 4214.560
## Ag2     4131 4214.379
## Ag3     4100 4189.446
# Dominancia
Dominancia <- dominance(merged_metagenomes, index = "all")
Dominancia
##            dbp        dmn absolute   relative     simpson core_abundance
## Ag1 0.02651136 0.04183593    40397 0.02651136 0.002837114      0.2763476
## Ag2 0.02891843 0.05012771    49917 0.02891843 0.003259692      0.2906512
## Ag3 0.02640574 0.04186274    39700 0.02640574 0.002898346      0.2796986
##          gini
## Ag1 0.8010570
## Ag2 0.8116390
## Ag3 0.8035928
# Rareza
Rareza <- rarity(merged_metagenomes, index = "all")
Rareza
##     log_modulo_skewness low_abundance rare_abundance
## Ag1            2.054063     0.7161781      0.7006737
## Ag2            2.052640     0.7008865      0.6776919
## Ag3            2.053428     0.7180931      0.6992699
# Cobertura
Cobertura <- coverage(merged_metagenomes, threshold = 0.5)
Cobertura
## Ag1 Ag2 Ag3 
## 231 212 228
# Desigualdad
Desigualdad <- inequality(merged_metagenomes)
Desigualdad
##       Ag1       Ag2       Ag3 
## 0.8010570 0.8116390 0.8035928
# Uniformidad
Uniformidad <- evenness(merged_metagenomes, "all")
Uniformidad
##       camargo    pielou    simpson      evar     bulla
## Ag1 0.2035280 0.8271414 0.08550968 0.1350859 0.3609202
## Ag2 0.1921895 0.8180177 0.07426228 0.1312996 0.3506728
## Ag3 0.2018350 0.8258311 0.08415229 0.1360249 0.3592546

Graficamos los resultados y calculamos la significancia estadística. Para esto usamos el paquete ggpubr que genera “publication-ready plots”, algo que siempre es deseable (ejecuta library(ggpubr)).

library("ggpubr")
# Generamos un objeto `phyloseq` sin taxa´s que tengan 0 reads
merged_metagenomes <- prune_taxa(taxa_sums(merged_metagenomes) > 0, merged_metagenomes)
# Calculamos los índices de diversidad
Alpha2 <- alpha(merged_metagenomes, index = "all")
# Y finalmente visualizamos la tabla de resultados
head(Alpha2)
##     observed    chao1 diversity_inverse_simpson diversity_gini_simpson
## Ag1     4122 4214.560                  352.4709              0.9971629
## Ag2     4131 4214.379                  306.7775              0.9967403
## Ag3     4100 4189.446                  345.0244              0.9971017
##     diversity_shannon diversity_fisher diversity_coverage evenness_camargo
## Ag1          6.885202         515.8116                231        0.2035280
## Ag2          6.811040         508.0487                212        0.1921895
## Ag3          6.869876         513.6507                228        0.2018350
##     evenness_pielou evenness_simpson evenness_evar evenness_bulla dominance_dbp
## Ag1       0.8271414       0.08550968     0.1350859      0.3609202    0.02651136
## Ag2       0.8180177       0.07426228     0.1312996      0.3506728    0.02891843
## Ag3       0.8258311       0.08415229     0.1360249      0.3592546    0.02640574
##     dominance_dmn dominance_absolute dominance_relative dominance_simpson
## Ag1    0.04183593              40397         0.02651136       0.002837114
## Ag2    0.05012771              49917         0.02891843       0.003259692
## Ag3    0.04186274              39700         0.02640574       0.002898346
##     dominance_core_abundance dominance_gini rarity_log_modulo_skewness
## Ag1                0.2763476      0.8010570                   2.054063
## Ag2                0.2906512      0.8116390                   2.052640
## Ag3                0.2796986      0.8035928                   2.053428
##     rarity_low_abundance rarity_rare_abundance
## Ag1            0.7161781             0.7006737
## Ag2            0.7008865             0.6776919
## Ag3            0.7180931             0.6992699

Metadatos con diversidad

Extraemos el metadata del objeto phyloseq.

merged_metagenomes.meta <- meta(merged_metagenomes)
head(merged_metagenomes.meta)
##     SAMPLE SITE  pH N_TOTAL P_OLSEN ORGANIC SITE_ID
## Ag1    Ag1    2 6.2    0.36    29.6    12.0     Ag1
## Ag2    Ag2    3 6.7    0.22    13.5    13.5     Ag2
## Ag3    Ag3    1 6.1    0.29     4.0    24.2     Ag3

Le agregamos la tabla de diversidad a la metadata

names(Alpha2)
##  [1] "observed"                   "chao1"                     
##  [3] "diversity_inverse_simpson"  "diversity_gini_simpson"    
##  [5] "diversity_shannon"          "diversity_fisher"          
##  [7] "diversity_coverage"         "evenness_camargo"          
##  [9] "evenness_pielou"            "evenness_simpson"          
## [11] "evenness_evar"              "evenness_bulla"            
## [13] "dominance_dbp"              "dominance_dmn"             
## [15] "dominance_absolute"         "dominance_relative"        
## [17] "dominance_simpson"          "dominance_core_abundance"  
## [19] "dominance_gini"             "rarity_log_modulo_skewness"
## [21] "rarity_low_abundance"       "rarity_rare_abundance"
merged_metagenomes.meta$Shannon <- Alpha2$diversity_shannon 
merged_metagenomes.meta$InverseSimpson <- Alpha2$diversity_inverse_simpson
merged_metagenomes.meta$Chao1 <- Alpha2$chao1
merged_metagenomes.meta$Observed <- Alpha2$observed
merged_metagenomes.meta$Fisher <- Alpha2$diversity_fisher

En este ejercicio nos interesa comparar la diversidad entre áreas. Recordemos que tenemos datos para tres áreas: Ag1, Ag2, Ag3. Necesitamos crear una lista de comparasiones de a pares para poder visualizar y calcular significancia estadística de manera simultánea.

# Obtenemos las variables desde nuestro objeto `phyloseq`
Area <- levels(as.factor(merged_metagenomes.meta$SITE_ID))
# Creamos una lista de lo que queremos comparar
ParesArea <- combn(seq_along(Area), 2, simplify = FALSE, FUN = function(i) Area[i])
ParesArea
## [[1]]
## [1] "Ag1" "Ag2"
## 
## [[2]]
## [1] "Ag1" "Ag3"
## 
## [[3]]
## [1] "Ag2" "Ag3"

Con la función ggviolin podemos generar un gráfico de violín en un solo paso de la siguiente forma.

# Cargamos una variable de color
# COLOR <- c("#a6cee3", "#b2df8a", "#fdbf6f")
Indices <- c("Observed","Shannon","Fisher","Chao1","InverseSimpson")
# Ciclo
for (i in Indices){
  # Prueba 
  # i <- "Observed"
  p <- ggviolin(merged_metagenomes.meta, x = "SITE_ID", y = i, 
                add = "boxplot", fill = "SITE_ID") + #, palette = COLOR)
    scale_fill_brewer(palette = "Set1") +
    # Evaluamos la significancia estadística
    stat_compare_means(comparisons = ParesArea)
  print(p)
}
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`
## Caused by error in `$<-.data.frame`:
## ! replacement has 1 row, data has 0
Gráfico de Violín de los diferentes índices de diversidad alfa $lpha$ y su significancia estadística.

Gráfico de Violín de los diferentes índices de diversidad alfa \(lpha\) y su significancia estadística.

## Warning: Groups with fewer than two data points have been dropped.
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`
## Caused by error in `$<-.data.frame`:
## ! replacement has 1 row, data has 0
Gráfico de Violín de los diferentes índices de diversidad alfa $lpha$ y su significancia estadística.

Gráfico de Violín de los diferentes índices de diversidad alfa \(lpha\) y su significancia estadística.

## Warning: Groups with fewer than two data points have been dropped.
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`
## Caused by error in `$<-.data.frame`:
## ! replacement has 1 row, data has 0
Gráfico de Violín de los diferentes índices de diversidad alfa $lpha$ y su significancia estadística.

Gráfico de Violín de los diferentes índices de diversidad alfa \(lpha\) y su significancia estadística.

## Warning: Groups with fewer than two data points have been dropped.
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`
## Caused by error in `$<-.data.frame`:
## ! replacement has 1 row, data has 0
Gráfico de Violín de los diferentes índices de diversidad alfa $lpha$ y su significancia estadística.

Gráfico de Violín de los diferentes índices de diversidad alfa \(lpha\) y su significancia estadística.

## Warning: Groups with fewer than two data points have been dropped.
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Warning in max(data$density): no non-missing arguments to max; returning -Inf
## Warning: Computation failed in `stat_ydensity()`
## Caused by error in `$<-.data.frame`:
## ! replacement has 1 row, data has 0
Gráfico de Violín de los diferentes índices de diversidad alfa $lpha$ y su significancia estadística.

Gráfico de Violín de los diferentes índices de diversidad alfa \(lpha\) y su significancia estadística.

Beta Diversidad

La diversidad beta \(\beta\) es la relación entre la diversidad de especies regional y local. La diversidad beta como una medida del recambio de especies enfatiza demasiado el papel de las especies raras, ya que la diferencia en la composición de especies entre dos sitios o comunidades probablemente refleje la presencia o ausencia de algunas especies raras en los ensamblajes. Existen diferentes medidas de similitud (o disimilitud, i.e., 1 - similitud) disponibles que nos permiten entender las relaciones entre nuestras muestras. En general todas producen matrices de distancia comparables. El paquete phyloseq ofrece un gran número de medidas de distancia. Las más populares son UniFrac y Weighted UniFrac (medidas que consideran filogenia) y otras independientes de filogenia como: Jaccard, Manhattan, Euclidian, Bray-Curtis, Canberra, etc. Por otra parte, la matriz de distancia resultante no se usa en aislación sino que en conjunto con algún método de ordinación o escalamiento multidimensional (ordination). De nuevo, phyloseq ofrece un gran número de métodos entre los cuales se encuentran: detrended y canonical correspondence analysis, Double Principal Coordinate Analysis, Non-metric MultiDimenstional Scaling, y MDS/PCoA. Probemos entonces hacer un análisis tipo PCoA con una matriz de distancia que considera las relaciones filogenéticas y otra que no.

# Con esto vemos todos los métodos (métricas de distancia) para calcular la beta diversidad
distanceMethodList
## $UniFrac
## [1] "unifrac"  "wunifrac"
## 
## $DPCoA
## [1] "dpcoa"
## 
## $JSD
## [1] "jsd"
## 
## $vegdist
##  [1] "manhattan"  "euclidean"  "canberra"   "bray"       "kulczynski"
##  [6] "jaccard"    "gower"      "altGower"   "morisita"   "horn"      
## [11] "mountford"  "raup"       "binomial"   "chao"       "cao"       
## 
## $betadiver
##  [1] "w"   "-1"  "c"   "wb"  "r"   "I"   "e"   "t"   "me"  "j"   "sor" "m"  
## [13] "-2"  "co"  "cc"  "g"   "-3"  "l"   "19"  "hk"  "rlb" "sim" "gl"  "z"  
## 
## $dist
## [1] "maximum"   "binary"    "minkowski"
## 
## $designdist
## [1] "ANY"
# Algunos métodos de ORDENACIÓN
Ordenacion <- c("NMDS", "MDS", "PCoA", "DPCoA")
# Algunas DISTANCIAS
Distancia <- c("bray", "unifrac", "dpcoa", "jsd", "jaccard")

#Para poder usar la variable año en la bandera "shape" la convertimos en factor
merged_metagenomes@sam_data@.Data[[6]] <- as.factor(merged_metagenomes@sam_data@.Data[[6]])
for (o in Ordenacion){
  # Prueba 
  # o <- "NMDS"
  for (d in Distancia){
    tryCatch({
      # Prueba 
      # d <- "bray"
      Ord <- ordinate(physeq = merged_metagenomes, method = o, distance = d)
      Beta <- plot_ordination(physeq = merged_metagenomes, ordination = Ord, 
                              color = "SITE_ID", shape = "YEAR") +
        theme_classic() +
        geom_point(size = 5) +
        geom_text(mapping = aes(label = SITE_ID), size = 4, vjust = 2, hjust = 1) + 
        geom_vline(xintercept = 0) +
        geom_hline(yintercept = 0) +
        labs(title = paste("Beta diversity with", o, "-", d, 
                           "for the Calakmul samples ShotGun")) +
               theme(legend.position = "right", 
                     legend.text = element_text(size=12, face="bold"),
                     legend.title = element_text(size = 12, face = "bold"),
                     axis.text = element_text(color = "black", size = 12, face = "bold"),
                     axis.title = element_text(color = "black", size = 12, face = "bold"),
                     plot.title = element_text(hjust = 0.5, size=14, face="bold")) +
               scale_color_brewer(palette = "Set1") + 
               scale_fill_brewer(palette = "Set1")
      print(Beta)
      pdf(paste0("Resultados/BetaDiversidad/PlotBeta", o, "-", d, ".pdf"), 
          width=10, height=7)
      print(Beta)
      dev.off()
      svg(paste0("Resultados/BetaDiversidad/PlotBeta", o, "-", d, ".svg"),
          width=10, height=7)
      print(Beta)
      dev.off()
    }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")}) 
  }
}
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0 
## Run 1 stress 0 
## ... Procrustes: rmse 0.02778263  max resid 0.02891375 
## Run 2 stress 0 
## ... Procrustes: rmse 0.3444346  max resid 0.4253598 
## Run 3 stress 0 
## ... Procrustes: rmse 0.06473031  max resid 0.07607818 
## Run 4 stress 0 
## ... Procrustes: rmse 0.02035538  max resid 0.02259239 
## Run 5 stress 0 
## ... Procrustes: rmse 0.06570887  max resid 0.07694945 
## Run 6 stress 0 
## ... Procrustes: rmse 0.1943052  max resid 0.2572746 
## Run 7 stress 0 
## ... Procrustes: rmse 0.1172204  max resid 0.1400096 
## Run 8 stress 0 
## ... Procrustes: rmse 0.05258035  max resid 0.0601127 
## Run 9 stress 0 
## ... Procrustes: rmse 0.3032284  max resid 0.4092741 
## Run 10 stress 0 
## ... Procrustes: rmse 0.07178567  max resid 0.08471621 
## Run 11 stress 0 
## ... Procrustes: rmse 0.3213424  max resid 0.4499814 
## Run 12 stress 0 
## ... Procrustes: rmse 0.3369762  max resid 0.47323 
## Run 13 stress 0 
## ... Procrustes: rmse 0.1224177  max resid 0.1528506 
## Run 14 stress 0 
## ... Procrustes: rmse 0.1830156  max resid 0.2293378 
## Run 15 stress 0 
## ... Procrustes: rmse 0.2391714  max resid 0.3252009 
## Run 16 stress 0 
## ... Procrustes: rmse 0.1224816  max resid 0.1484967 
## Run 17 stress 0 
## ... Procrustes: rmse 0.1958854  max resid 0.2594873 
## Run 18 stress 0 
## ... Procrustes: rmse 0.1352705  max resid 0.1625515 
## Run 19 stress 0 
## ... Procrustes: rmse 0.3467075  max resid 0.4864014 
## Run 20 stress 0 
## ... Procrustes: rmse 0.07094986  max resid 0.0839901 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress < smin
## Warning in metaMDS(veganifyOTU(physeq), distance, ...): stress is (nearly) zero:
## you may have insufficient data
## Warning in postMDS(out$points, dis, plot = max(0, plot - 1), ...): skipping
## half-change scaling: too few points below threshold
## Warning in plot_ordination(physeq = merged_metagenomes, ordination = Ord, :
## Shape variable was not found in the available data you provided.No shape mapped.
## ERROR : phy_tree slot is empty. 
## ERROR : phy_tree slot is empty. 
## Run 0 stress 0 
## Run 1 stress 0 
## ... Procrustes: rmse 0.2344238  max resid 0.2937118 
## Run 2 stress 0 
## ... Procrustes: rmse 0.1935433  max resid 0.2591064 
## Run 3 stress 0 
## ... Procrustes: rmse 0.3277101  max resid 0.3423276 
## Run 4 stress 0 
## ... Procrustes: rmse 0.2907152  max resid 0.3136829 
## Run 5 stress 0 
## ... Procrustes: rmse 0.2350237  max resid 0.2776408 
## Run 6 stress 0 
## ... Procrustes: rmse 0.2392581  max resid 0.2671299 
## Run 7 stress 0 
## ... Procrustes: rmse 0.1867016  max resid 0.2618725 
## Run 8 stress 0 
## ... Procrustes: rmse 0.2139299  max resid 0.2524363 
## Run 9 stress 0 
## ... Procrustes: rmse 0.1946998  max resid 0.2310102 
## Run 10 stress 0 
## ... Procrustes: rmse 0.2921048  max resid 0.3224547 
## Run 11 stress 0 
## ... Procrustes: rmse 0.1562071  max resid 0.2136159 
## Run 12 stress 0 
## ... Procrustes: rmse 0.2650681  max resid 0.3206244 
## Run 13 stress 0 
## ... Procrustes: rmse 0.169596  max resid 0.2018941 
## Run 14 stress 0 
## ... Procrustes: rmse 0.2484225  max resid 0.2982738 
## Run 15 stress 0 
## ... Procrustes: rmse 0.04849595  max resid 0.06031097 
## Run 16 stress 0 
## ... Procrustes: rmse 0.1296971  max resid 0.1542298 
## Run 17 stress 0 
## ... Procrustes: rmse 0.2222494  max resid 0.253172 
## Run 18 stress 0 
## ... Procrustes: rmse 0.07582651  max resid 0.09313853 
## Run 19 stress 0 
## ... Procrustes: rmse 0.02437069  max resid 0.03032981 
## Run 20 stress 0 
## ... Procrustes: rmse 0.04372944  max resid 0.05389333 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress < smin
## Warning in metaMDS(ps.dist): stress is (nearly) zero: you may have insufficient
## data

## Warning in metaMDS(ps.dist): Shape variable was not found in the available data
## you provided.No shape mapped.
ïndices de Diversidad Beta $eta$.

ïndices de Diversidad Beta \(eta\).

## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0 
## Run 1 stress 0 
## ... Procrustes: rmse 0.2480075  max resid 0.3385209 
## Run 2 stress 0 
## ... Procrustes: rmse 0.1603973  max resid 0.2038058 
## Run 3 stress 0 
## ... Procrustes: rmse 0.2789778  max resid 0.3863598 
## Run 4 stress 0 
## ... Procrustes: rmse 0.1478214  max resid 0.18778 
## Run 5 stress 0 
## ... Procrustes: rmse 0.2720961  max resid 0.3705072 
## Run 6 stress 0 
## ... Procrustes: rmse 0.3456619  max resid 0.4254321 
## Run 7 stress 0 
## ... Procrustes: rmse 0.1029492  max resid 0.1256955 
## Run 8 stress 0 
## ... Procrustes: rmse 0.1500741  max resid 0.1910323 
## Run 9 stress 0 
## ... Procrustes: rmse 0.0219684  max resid 0.02314271 
## Run 10 stress 0 
## ... Procrustes: rmse 0.2387852  max resid 0.3208074 
## Run 11 stress 0 
## ... Procrustes: rmse 0.2281143  max resid 0.2732617 
## Run 12 stress 0 
## ... Procrustes: rmse 0.3162651  max resid 0.3952995 
## Run 13 stress 0 
## ... Procrustes: rmse 0.06210309  max resid 0.0624079 
## Run 14 stress 0 
## ... Procrustes: rmse 0.04787707  max resid 0.04892281 
## Run 15 stress 0 
## ... Procrustes: rmse 0.3288374  max resid 0.4626144 
## Run 16 stress 0 
## ... Procrustes: rmse 0.1194385  max resid 0.1475403 
## Run 17 stress 0 
## ... Procrustes: rmse 0.2794117  max resid 0.3515494 
## Run 18 stress 0 
## ... Procrustes: rmse 0.3443004  max resid 0.4315653 
## Run 19 stress 0 
## ... Procrustes: rmse 0.02595202  max resid 0.02774001 
## Run 20 stress 0 
## ... Procrustes: rmse 0.1491031  max resid 0.1895109 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress < smin
## Warning in metaMDS(veganifyOTU(physeq), distance, ...): stress is (nearly) zero:
## you may have insufficient data
## Warning in postMDS(out$points, dis, plot = max(0, plot - 1), ...): skipping
## half-change scaling: too few points below threshold
## Warning in plot_ordination(physeq = merged_metagenomes, ordination = Ord, :
## Shape variable was not found in the available data you provided.No shape mapped.
ïndices de Diversidad Beta $eta$.

ïndices de Diversidad Beta \(eta\).

## Warning in plot_ordination(physeq = merged_metagenomes, ordination = Ord, :
## Shape variable was not found in the available data you provided.No shape mapped.
ïndices de Diversidad Beta $eta$.

ïndices de Diversidad Beta \(eta\).

## ERROR : phy_tree slot is empty. 
## ERROR : phy_tree slot is empty.
## Warning in plot_ordination(physeq = merged_metagenomes, ordination = Ord, :
## Shape variable was not found in the available data you provided.No shape mapped.
ïndices de Diversidad Beta $eta$.

ïndices de Diversidad Beta \(eta\).

## Warning in plot_ordination(physeq = merged_metagenomes, ordination = Ord, :
## Shape variable was not found in the available data you provided.No shape mapped.
ïndices de Diversidad Beta $eta$.

ïndices de Diversidad Beta \(eta\).

## Warning in plot_ordination(physeq = merged_metagenomes, ordination = Ord, :
## Shape variable was not found in the available data you provided.No shape mapped.
ïndices de Diversidad Beta $eta$.

ïndices de Diversidad Beta \(eta\).

## ERROR : phy_tree slot is empty. 
## ERROR : phy_tree slot is empty.
## Warning in plot_ordination(physeq = merged_metagenomes, ordination = Ord, :
## Shape variable was not found in the available data you provided.No shape mapped.
ïndices de Diversidad Beta $eta$.

ïndices de Diversidad Beta \(eta\).

## Warning in plot_ordination(physeq = merged_metagenomes, ordination = Ord, :
## Shape variable was not found in the available data you provided.No shape mapped.
ïndices de Diversidad Beta $eta$.

ïndices de Diversidad Beta \(eta\).

ïndices de Diversidad Beta $eta$.

ïndices de Diversidad Beta \(eta\).

## ERROR : phy_tree slot is empty. 
## ERROR : phy_tree slot is empty. 
## ERROR : phy_tree slot is empty. 
## ERROR : phy_tree slot is empty. 
## ERROR : phy_tree slot is empty.

Nota que los gráficos de dispersión donde se visualiza este tipo de análisis están escalados según la cantidad de variación que los ejes explican. En general lo que estos métodos pretenden hacer es tratar de encontrar el menor número de vectores matemáticos que maximicen la separación entre las muestras (puntos en el gráfico). Esto nace de la imposibilidad de graficar eficientemente datos multidimencionales. Volviendo a los ejes, estos no suman 100% porque hay otros ejes que no estamos usando para graficar y que contribuyen con el resto de la variación. Al graficarlos de manera simétrica distorsionaríamos las relaciones entre los puntos, especialmente si estamos comparando dos o más gráficos. En específico para comunidades microbianas, el método de doble análisis de coordenadas principales o (DPCoA) es muy apropiado porque analiza conjuntamente dos tipos de datos: una tabla de disimilitud que representa diferencias entre especies y una matriz de abundancia que representa la distribución de especies entre las comunidades. El resultado final es un ensamble del espacio multidimensional que correlacionan las especies con las comunidades. Ahora exploremos escalamiento multidimensional usando un método reciente conocido como t-SNE o t-Distributed Stochastic Neighbor Embedding. t-SNE difiere de otros métodos en que hace énfasis en las distancias locales en vez de las distancias globales, de esta forma generando una mayor resolución o separación entre los puntos o muestras.

Distancia <- c("bray", "unifrac", "dpcoa", "jsd")
for (d in Distancia){
  tryCatch({
    # Prueba 
    # d <- "bray"
    tSNE <- tsne_phyloseq(physeq = merged_metagenomes, distance= d, perplexity = 8, 
                          verbose=0, rng_seed = 3901)
    beta <- plot_tsne_phyloseq(physeq = merged_metagenomes, tSNE, color = "SITE_ID") +
      theme_classic() +
      geom_point(size = 5) +
      geom_text(mapping = aes(label = SITE_ID), size = 4, vjust = 2, hjust = 1) + 
      geom_vline(xintercept = 0) +
      geom_hline(yintercept = 0) +
      labs(title = paste("Beta diversity with tSNE", "-", d, "for the Calakmul samples ShotGun")) +
      theme(legend.position = "right", legend.text = element_text(size=12, face="bold"),
            legend.title = element_text(size = 12, face = "bold"),
            axis.text = element_text(color = "black", size = 12, face = "bold"),
            axis.title = element_text(color = "black", size = 12, face = "bold"),
            plot.title = element_text(hjust = 0.5, size=14, face="bold")) +
      scale_color_brewer(palette = "Set1") + 
      scale_fill_brewer(palette = "Set1")
    print(beta)
    pdf(paste0("Resultados/BetaDiversidad/PlotBetatSNE-", d, ".pdf"), 
        width=10, height=7)
    print(beta)
    dev.off()
    svg(paste0("Resultados/BetaDiversidad/PlotBetatSNE-", d, ".svg"),
        width=10, height=7)
    print(beta)
  }, error=function(e){cat("ERROR :",conditionMessage(e), "\n")}) 
}
## ERROR : Failed to tune perplexity at observation 1. Try using a smaller perlexity. 
## ERROR : Phylogenetic tree not available in physeq obejct. distance = unifrac does not work. 
## ERROR : phy_tree slot is empty. 
## ERROR : Failed to tune perplexity at observation 1. Try using a smaller perlexity.

Análisis de abundancias y visualizaciones

Una manera muy eficiente de obtener una visión general de la composición de las muestras es a través de un gráfico de barras apiladas.

# Necesitamos obtener las taxa más abundantes, en este caso el top 10
top10 <- get_top_taxa(physeq_obj = merged_metagenomes, n = 10, relative = T,
                       discard_other = T, other_label = "Other")
## Warning: 'get_top_taxa' is deprecated.
## Use 'top_taxa' instead.
## Use 'nested_top_taxa' instead.
## See help("Deprecated")
# Ya que no todas las taxa fueron clasificadas a nivel de especie, generamos etiquetas compuestas de distintos rangos taxonómicos para el gráfico
top10 <- name_taxa(top10, label = "", species = F, other_label = "Other")
## Warning: 'name_taxa' is deprecated.
## Use 'name_na_taxa' instead.
## Use 'label_duplicate_taxa' instead.
## See help("Deprecated")
# Finalmente graficamos
ptop10 <- fantaxtic_bar(top10, color_by = "Family", label_by = "Genus", 
                        facet_by = NULL, grid_by = NULL, other_color = "Grey") +
  theme_classic() +
  labs(title="Top 10 genus present in Calakmul samples obtained by ShotGun")+
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text.x = element_text(angle = 45, vjust = 1.0, hjust=1),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold")) + 
  scale_fill_brewer(palette = "Set3") 
## Warning: 'fantaxtic_bar' is deprecated.
## Use 'fantaxtic' instead.
## See help("Deprecated")
## Warning: 'gen_palette' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")
## Warning: 'gen_colors' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")
## Warning: 'shuffle_colors' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")
## Warning: 'gen_shades_tints' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")

## Warning: 'gen_shades_tints' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")

## Warning: 'gen_shades_tints' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")

## Warning: 'gen_shades_tints' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")

## Warning: 'gen_shades_tints' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")

## Warning: 'gen_shades_tints' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")

## Warning: 'gen_shades_tints' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")
## Warning: 'gen_uniq_lbls' is deprecated.
## Use 'name_na_taxa' instead.
## Use 'label_duplicate_taxa' instead.
## See help("Deprecated")
ptop10
Top 10 de géneros más abundantes y clasificado por *Familia*.

Top 10 de géneros más abundantes y clasificado por Familia.

pdf("Resultados/PlotTop10genus.pdf", width=10, height=7)
ptop10
dev.off()
## png 
##   2
svg("Resultados/PlotTop10genus.svg", width=10, height=7)
ptop10
dev.off()
## png 
##   2
# facet_by
ptop10.1 <- fantaxtic_bar(top10, color_by = "Family", label_by = "Genus", 
                        facet_by = "SITE_ID", grid_by = NULL, other_color = "Grey") +
  theme_classic() +
  labs(title="Top 10 genus present in Calakmul samples obtained by ShotGun")+
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text.x = element_text(angle = 45, vjust = 1.0, hjust=1),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold")) + 
  scale_fill_brewer(palette = "Set3") 
## Warning: 'fantaxtic_bar' is deprecated.
## Use 'fantaxtic' instead.
## See help("Deprecated")
## Warning: 'gen_palette' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")
## Warning: 'gen_colors' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")
## Warning: 'shuffle_colors' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")
## Warning: 'gen_shades_tints' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")

## Warning: 'gen_shades_tints' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")

## Warning: 'gen_shades_tints' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")

## Warning: 'gen_shades_tints' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")

## Warning: 'gen_shades_tints' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")

## Warning: 'gen_shades_tints' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")

## Warning: 'gen_shades_tints' is deprecated.
## Use 'nested_palette' instead.
## See help("Deprecated")
## Warning: 'gen_uniq_lbls' is deprecated.
## Use 'name_na_taxa' instead.
## Use 'label_duplicate_taxa' instead.
## See help("Deprecated")
ptop10.1
Top 10 de géneros más abundantes y clasificado por *Familia*.

Top 10 de géneros más abundantes y clasificado por Familia.

pdf("Resultados/PlotTop10genus1.pdf", width=10, height=7)
ptop10.1
dev.off()
## png 
##   2
svg("Resultados/PlotTop10genus1.svg", width=10, height=7)
ptop10.1
dev.off()
## png 
##   2

Diferentes visualizaciones de abundancias

Veamos ahora otras herramientas que nos permiten examinar la composición de estas comunidades microbianas. El paquete ampvis2, desarrollado por Mads Albertsen y Kasper Skytte Andersen, nos permite hacer precisamente esto. Primero transformemos el objeto phyloseq con el cual hemos estado trabajando en un objeto ampvis2.

library(ampvis2)
# Necesitamos extraer la tabla de read counts y la tabla de taxonomía del objeto merged_metagenomes
# Generamos una copia para no sobreescribir merged_metagenomes
obj <- merged_metagenomes
# Cambiamos la orientación de la otu_table
otu_table(obj) <- t(otu_table(obj))
# Extraemos las tablas
otutable <- data.frame(OTU = rownames(t(phyloseq::otu_table(obj)@.Data)),
                       t(phyloseq::otu_table(obj)@.Data),
                       phyloseq::tax_table(obj)@.Data,
                       check.names = FALSE)
# Extraemos la metadada
metadata <- data.frame(phyloseq::sample_data(obj), 
                       check.names = FALSE)

Ampvis2 requiere - Los rangos taxonómicos sean y vayan de Kingdom a Species - La primera columna del metadata sea el identificador de cada muestra

# Entonces duplicamos la columna Género y le cambiamos el nombre a Especie
otutable$Species = otutable$Genus

# Finalmente generamos el objeto ampvis
AV2 <- amp_load(otutable, metadata)

Ahora echemos un vistazo a la estructura de las comunidades utilizando “rank abundance curves”.

# No logaritmico
AV2.abun <- amp_rankabundance(AV2, group_by = "SITE_ID", showSD = T, log10_x = F) +
  theme_classic() +
  labs(title="Relative abundance by number of accumulated readings")+
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold")) +
  geom_vline(xintercept = 1000) +
  geom_hline(yintercept = 90) +
  geom_vline(xintercept = 25) +
  geom_hline(yintercept = 10) +
  scale_fill_brewer(palette = "Set1") +
  scale_color_brewer(palette = "Set1")
AV2.abun
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
Estructura de comunidades por "Curvas de abundancia por rango, logarítmico y no logarítimico.

Estructura de comunidades por “Curvas de abundancia por rango, logarítmico y no logarítimico.

pdf("Resultados/AbundanciasAV2/PlotAmpAbun.pdf", width=10, height=7)
AV2.abun
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
dev.off()
## png 
##   2
svg("Resultados/AbundanciasAV2/PlotAmpAbun.svg", width=10, height=7)
AV2.abun
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
dev.off()
## png 
##   2
# Logaritmico
AV2.abunlog <- amp_rankabundance(AV2, group_by = "SITE_ID", showSD = T, log10_x = T) +
  theme_classic() +
  labs(title="Relative abundance by number of accumulated readings")+
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold")) +
  geom_vline(xintercept = 1000) +
  geom_hline(yintercept = 90) +
  geom_vline(xintercept = 25) +
  geom_hline(yintercept = 10) +
  scale_fill_brewer(palette = "Set1") +
  scale_color_brewer(palette = "Set1")
AV2.abunlog
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
Estructura de comunidades por "Curvas de abundancia por rango, logarítmico y no logarítimico.

Estructura de comunidades por “Curvas de abundancia por rango, logarítmico y no logarítimico.

pdf("Resultados/AbundanciasAV2/PlotAmpAbunLog.pdf", width=10, height=7)
AV2.abunlog
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
dev.off()
## png 
##   2
svg("Resultados/AbundanciasAV2/PlotAmpAbunLog.svg", width=10, height=7)
AV2.abunlog
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
dev.off()
## png 
##   2

El gráfico nos muestra que en la medida que vamos sumando las taxa de mayor a menor abundancia (Rank Abundance) la abundancia de reads cumulativa va aumentando. Lo importante de observar es la forma de la curva. Una curva que sube rápidamente nos indicaría que las comunidades están dominadas por unas cuantas taxa. El objeto AV2 tiene un total de 2220 OTUs, los que están entre un rango de 80 y 1000 están tomando el 90 % de las lecturas,lo que sugiere que el total de las lecturas está determinado por la mitad de todos los OTUs. Ahora veremos cuales son mediante un mapa de calor

# Phylum y Género
AmpHM_GP <- amp_heatmap(AV2, 
            group_by = "SITE_ID", facet_by = "SITE_ID", plot_values = TRUE,
            tax_show = 15, tax_aggregate = "Genus", tax_add = "Phylum",
            # color_vector = c("lightblue3", "olivedrab3"),
            plot_colorscale = "sqrt", plot_legendbreaks = c(0, 0.5, 1, 1.5, 2, 2.5)) +
  labs(title="Relative abundance for Phylum and genus")+
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold"))
AmpHM_GP
Mapas de calor a nivel de "*Phylum*" y/o "*Genus*"

Mapas de calor a nivel de “Phylum” y/o “Genus

pdf("Resultados/AbundanciasAV2/PlotAmpHM_GP.pdf", width=10, height=7)
AmpHM_GP 
dev.off()
## png 
##   2
svg("Resultados/AbundanciasAV2/PlotAmpHM_GP.svg", width=10, height=7)
AmpHM_GP 
dev.off()
## png 
##   2
# Phylum
AmpHM_P <- amp_heatmap(AV2, 
            group_by = "SITE_ID", facet_by = "SITE_ID", plot_values = TRUE,
            tax_show = 15, tax_aggregate = "Phylum",
            # color_vector = c("lightblue3", "olivedrab3"),
            plot_colorscale = "sqrt", plot_legendbreaks = c(0, 1, 5, 10, 20, 35)) +
  labs(title="Relative abundance for Phylum")+
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold"))
AmpHM_P
Mapas de calor a nivel de "*Phylum*" y/o "*Genus*"

Mapas de calor a nivel de “Phylum” y/o “Genus

pdf("Resultados/AbundanciasAV2/PlotAmpHM_P.pdf", width=10, height=7)
AmpHM_P 
dev.off()
## png 
##   2
svg("Resultados/AbundanciasAV2/PlotAmpHM_P.svg", width=10, height=7)
AmpHM_P 
dev.off()
## png 
##   2
# Género
AmpHM_G <- amp_heatmap(AV2, 
            group_by = "SITE_ID", facet_by = "SITE_ID", plot_values = TRUE,
            tax_show = 15, tax_aggregate = "Genus",
            # color_vector = c("lightblue3", "olivedrab3"),
            plot_colorscale = "sqrt", plot_legendbreaks = c(0, 1, 5, 10, 20, 35)) +
  labs(title="Relative abundance for Genus")+
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold"))
AmpHM_G
Mapas de calor a nivel de "*Phylum*" y/o "*Genus*"

Mapas de calor a nivel de “Phylum” y/o “Genus

pdf("Resultados/AbundanciasAV2/PlotAmpHM_G.pdf", width=10, height=7)
AmpHM_G 
dev.off()
## png 
##   2
svg("Resultados/AbundanciasAV2/PlotAmpHM_G.svg", width=10, height=7)
AmpHM_G
dev.off()
## png 
##   2

A nivel de “phylum” el más abundante en los tres sitios es Actinobacteriota con porcentajes de 27 a 36 % seguido de Proteobacteria con porcentajes alrededor de 18 - 25 % y en tercer lugar Acidobacteria con porcentajes que van de 14 a 19 %. Respecto al género, son diferentes géneros quienes dominan en cada uno de los sitios, Solirubrobacter es más abundante en Ag1, Haliangium en Ag2 y Nocardioides en Ag3. También podemos realizar una visualización similar pero usando Box Plots.

# Phylum y Género
AmpBP_PG <- amp_boxplot(AV2, group_by = "SITE_ID", tax_show = 15, 
                        tax_aggregate = "Genus", tax_add = "Phylum",
                        adjust_zero = T, plot_log = T) +
  labs(title="Relative abundance for Phylum and genus")+
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold")) +
  scale_fill_brewer(palette = "Set1") +
  scale_color_brewer(palette = "Set1")
AmpBP_PG
Boxplot de calor a nivel de "*Phylum*" y/o "*Genus*"

Boxplot de calor a nivel de “Phylum” y/o “Genus

pdf("Resultados/AbundanciasAV2/PlotAmpBP_PG.pdf", width=10, height=7)
AmpBP_PG
dev.off()
## png 
##   2
svg("Resultados/AbundanciasAV2/PlotAmpBP_PG.svg", width=10, height=7)
AmpBP_PG
dev.off()
## png 
##   2
# Phylum 
AmpBP_P <- amp_boxplot(AV2, group_by = "SITE_ID", tax_show = 15, 
                        tax_aggregate = "Phylum", adjust_zero = T, plot_log = T) +
  labs(title="Relative abundance for Phylum")+
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold")) +
  scale_fill_brewer(palette = "Set1") +
  scale_color_brewer(palette = "Set1")
AmpBP_P
Boxplot de calor a nivel de "*Phylum*" y/o "*Genus*"

Boxplot de calor a nivel de “Phylum” y/o “Genus

pdf("Resultados/AbundanciasAV2/PlotAmpBP_P.pdf", width=10, height=7)
AmpBP_P
dev.off()
## png 
##   2
svg("Resultados/AbundanciasAV2/PlotAmpBP_P.svg", width=10, height=7)
AmpBP_P
dev.off()
## png 
##   2
# Género
AmpBP_G <- amp_boxplot(AV2, group_by = "SITE_ID", tax_show = 15, 
                        tax_aggregate = "Genus", adjust_zero = T, plot_log = T) +
  labs(title="Relative abundance for Genus")+
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank(),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        plot.title = element_text(hjust = 0.5, size=14, face="bold")) +
  scale_fill_brewer(palette = "Set1") +
  scale_color_brewer(palette = "Set1")
AmpBP_G
Boxplot de calor a nivel de "*Phylum*" y/o "*Genus*"

Boxplot de calor a nivel de “Phylum” y/o “Genus

pdf("Resultados/AbundanciasAV2/PlotAmpBP_G.pdf", width=10, height=7)
AmpBP_G
dev.off()
## png 
##   2
svg("Resultados/AbundanciasAV2/PlotAmpBP_G.svg", width=10, height=7)
AmpBP_G
dev.off()
## png 
##   2

Veamos ahora si es que algunos de estos microorganismos están compartidos entre todas las muestras. Para esto debemos calcular el core microbiome o el conjunto de taxa compartidas entre un cierto umbral porcentual de muestras y de prevalencia intra-muestra.

AmpCore <- amp_core(AV2, group_by = "SITE_ID", core_pct = 90, margin_plots = "xy",
         margin_plot_values_size = 3, widths = c(3, 1), heights = c(1, 3))
AmpCore
Plotcore a nivel de "*Phylum*" y/o "*Genus*"

Plotcore a nivel de “Phylum” y/o “Genus

pdf("Resultados/AbundanciasAV2/PlotAmpCore.pdf", width=10, height=7)
AmpCore
dev.off()
## png 
##   2
svg("Resultados/AbundanciasAV2/PlotAmpCore.svg", width=10, height=7)
AmpCore
dev.off()
## png 
##   2

Y visto de otra manera en un diagrama de Venn.

# cut_a % OTUs que debajo de esta abundancia se excluyen
# cut_f % OTUs que superior a este es la frecuencia que consideran el núcleo
AmpVenn <- amp_venn(AV2, group_by = "SITE_ID", cut_a = 0.10, cut_f = 20, text_size = 4)
AmpVenn
Diagrama de Venn del Core taxonómico de cada zona

Diagrama de Venn del Core taxonómico de cada zona

pdf("Resultados/AbundanciasAV2/PlotAmpVenn.pdf", width=10, height=7)
AmpVenn
dev.off()
## png 
##   2
svg("Resultados/AbundanciasAV2/PlotAmpVenn.svg", width=10, height=7)
AmpVenn
dev.off()
## png 
##   2

Pheatmap

Para obtener los taxones dominantes, solo mantendremos las OTU cuyas lecturas sean mayor a 100 de manera abosluta o que estén presentes en un porcentaje mayor a 1 % de lecturas relativas

Es necesario que nuestro objeto phyloseq contenga la tabla de “otu_table” la siguiente estructura - Filas con OTUs - Columnas con Muestras

head(merged_metagenomes@otu_table@.Data) # Ya está como la queremos
##           Ag1   Ag2   Ag3
## 374     23351 25994 23239
## 1437360 15340 16278 15959
## 1274631 13148 14339 12890
## 2057741  6885 10364  6922
## 722472   5475  5579  5395
## 1355477  5236  8097  5099
summary(merged_metagenomes@otu_table@.Data) # Resumen estadístico
##       Ag1               Ag2               Ag3         
##  Min.   :    0.0   Min.   :    0.0   Min.   :    0.0  
##  1st Qu.:    9.0   1st Qu.:    9.0   1st Qu.:    8.0  
##  Median :   41.0   Median :   40.0   Median :   38.0  
##  Mean   :  361.3   Mean   :  409.3   Mean   :  356.5  
##  3rd Qu.:  314.0   3rd Qu.:  339.0   3rd Qu.:  309.0  
##  Max.   :40397.0   Max.   :49917.0   Max.   :39700.0

Debido a que el análisis a nivel de “Specie” no es tan confiable, nos quedamos a nivel de género y además esto hace que los del mismo género los una en uno solo

# Cortamos los OTUs a nivel de GENUS
merged_metagenomes2 <- tax_glom(merged_metagenomes,taxrank = rank_names(merged_metagenomes)[6]) 
summary(merged_metagenomes2@otu_table@.Data)
##       Ag1                Ag2                Ag3        
##  Min.   :     0.0   Min.   :     0.0   Min.   :     0  
##  1st Qu.:    28.0   1st Qu.:    28.0   1st Qu.:    27  
##  Median :   149.5   Median :   141.0   Median :   138  
##  Mean   :  1400.5   Mean   :  1586.5   Mean   :  1382  
##  3rd Qu.:   866.8   3rd Qu.:   876.5   3rd Qu.:   851  
##  Max.   :149834.0   Max.   :191644.0   Max.   :149268

De aquí en adelante, trabajamos con valores Relativos_S y con valores Absolutos_S

# Absolutos_S
# Filtramos aquellos géneros que tienen una media mínima de 100 lecturas absolutas
psd.abs.100 <- filter_taxa(merged_metagenomes2, function(x) mean(x) > 10000.0, TRUE)

# Relativos_S
psd.rel <- transform_sample_counts(merged_metagenomes2, function(x) x*100 / sum(x))
# Filtramos aquellos géneros que tiene al menos un 1 % de lecturas relativas 
psd.rel.1 <- filter_taxa(psd.rel, function(x) mean(x) > 1.0, TRUE)

El paquete pheatmap necesista un dataframe como entrada, asi que haremos unos desde el objeto phyloseq

#Absolutos_S
psd.abs.100
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 27 taxa and 3 samples ]
## sample_data() Sample Data:       [ 3 samples by 7 sample variables ]
## tax_table()   Taxonomy Table:    [ 27 taxa by 7 taxonomic ranks ]
psd.abs.df <- as.data.frame(psd.abs.100@otu_table@.Data)
dim(psd.abs.df)
## [1] 27  3
rownames(psd.abs.df) <- make.names(psd.abs.100@tax_table@.Data[,6], unique = T)
str(psd.abs.df)
## 'data.frame':    27 obs. of  3 variables:
##  $ Ag1: num  103965 18866 16562 21430 19431 ...
##  $ Ag2: num  127690 19678 25139 23185 25729 ...
##  $ Ag3: num  104646 19052 16747 21702 20759 ...
#Relativos_S
psd.rel.1
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 18 taxa and 3 samples ]
## sample_data() Sample Data:       [ 3 samples by 7 sample variables ]
## tax_table()   Taxonomy Table:    [ 18 taxa by 7 taxonomic ranks ]
psd.rel.df <- as.data.frame(psd.rel.1@otu_table@.Data)
dim(psd.rel.df)
## [1] 18  3
rownames(psd.rel.df) <- make.names(psd.rel.1@tax_table@.Data[,6], unique = T)
str(psd.rel.df)
## 'data.frame':    18 obs. of  3 variables:
##  $ Ag1: num  6.82 1.24 1.09 1.41 1.28 ...
##  $ Ag2: num  7.4 1.14 1.46 1.34 1.49 ...
##  $ Ag3: num  6.96 1.27 1.11 1.44 1.38 ...

Ahora, necesitamos un nuevo marco de datos donde asignaremos los metadatos que queremos visualizar en el pheatmap

meta <- data.frame(Site = metadata$SITE_ID, row.names = rownames(metadata))
#meta$ID <- ID=metadata$ID
#meta$pH <- metadata$pH
#meta$N_Total <- metadata$N_TOTAL
#meta$P_Olsen <- metadata$P_OLSEN
#meta$Organic <- metadata$ORGANIC
# meta$Site <- metadata$SITE_ID
str(meta)
## 'data.frame':    3 obs. of  1 variable:
##  $ Site: chr  "Ag1" "Ag2" "Ag3"

Tambien generaremos una lista de colores

# Colores del metadatos
Color <- list(#ID=c("Ag2_Sep17"="antiquewhite", "Ag3_Sep17"="cadetblue1",
                  #    "Ag1_Sep17"="chartreuse", "Ag3_Sep18"="azure4", 
                  #    "Ag3_Jun19"="coral", "Ag2_Sep19"="bisque3",
                  #    "Ag2_Sep18"="bisque3", "Ag3_Sep19"="darkcyan",
                  #    "Ag1_Sep18"="darkmagenta", "Ag1_Sep19"="deepskyblue3"),
                 # Year = c("2017"="deepskyblue3", "2018"="indianred4", "2019"="purple4"),
                 #pH=c("6.7"="#bababa", "6.1"="brown4", "6.2"="#377eb8"),
                 #N_Total = c("0.22"="darkolivegreen", "0.29"="#984ea3", "0.36"="darkgoldenrod3"),
                 #P_Oolse = c("13.5"="darkred", "4.0"="deeppink4", "29.6"="aquamarine4"),
                 #Organic = c("13.5"="darkseagreen4", "24.2"="gold4", "12.0"="gray"),
                 Site = c("Ag1"="lightseagreen", "Ag2"="midnightblue", "Ag3"="mediumvioletred"))
# Color Pheatmap
ColorHPM <- c("#ffffff","#fff7fb","#ece2f0","#d0d1e6","#a6bddb","#67a9cf","#3690c0",
                    "#02818a","#016c59","#014636")

Corremos pheatmap

# Aboslutos
PHMap.Abs <- pheatmap(psd.abs.df, color = ColorHPM, cluster_cols = TRUE, 
                      cutree_cols = 3, cutree_rows = 4, border_color ="#000000",
                      annotation_col = meta, annotation_colors = Color)
PHMap.Abs
Mapas de calor usando Pheatmap con valores Absolutos_S y valores Relativos_S.

Mapas de calor usando Pheatmap con valores Absolutos_S y valores Relativos_S.

pdf("Resultados/PlotPHMap.Abs.pdf", width=20, height=10)
PHMap.Abs
dev.off()
## png 
##   2
svg("Resultados/PlotPHMap.Abs.svg", width=20, height=10)
PHMap.Abs
dev.off()
## png 
##   2
# Relativos_S
PHMap.Rel <- pheatmap(psd.rel.df, color = ColorHPM, cluster_cols = TRUE, 
                      cutree_cols = 3, cutree_rows = 4, border_color ="#000000",
                      annotation_col = meta, annotation_colors = Color)
PHMap.Rel
Mapas de calor usando Pheatmap con valores Absolutos_S y valores Relativos_S.

Mapas de calor usando Pheatmap con valores Absolutos_S y valores Relativos_S.

pdf("Resultados/PlotPHMap.Rel.pdf", width=20, height=10)
PHMap.Rel
dev.off()
## png 
##   2
svg("Resultados/PlotPHMap.Rel.svg", width=20, height=10)
PHMap.Rel
dev.off()
## png 
##   2

Otros gráficos

Vamos a crear un data.frame con el que trabajaremos

# Cortamos y agrupamos a nivel de género
Genero <- tax_glom(merged_metagenomes, taxrank = 'Genus')

# Absolutos_S
# Generamos el data.frame
V.Absolutos_S <- psmelt(Genero)
str(V.Absolutos_S)
## 'data.frame':    3264 obs. of  16 variables:
##  $ OTU      : chr  "1883" "1883" "1883" "374" ...
##  $ Sample   : chr  "Ag2" "Ag1" "Ag3" "Ag2" ...
##  $ Abundance: num  191644 149834 149268 127690 104646 ...
##  $ SAMPLE   : chr  "Ag2" "Ag1" "Ag3" "Ag2" ...
##  $ SITE     : chr  "3" "2" "1" "3" ...
##  $ pH       : num  6.7 6.2 6.1 6.7 6.1 6.2 6.1 6.7 6.2 6.2 ...
##  $ N_TOTAL  : num  0.22 0.36 0.29 0.22 0.29 0.36 0.29 0.22 0.36 0.36 ...
##  $ P_OLSEN  : num  13.5 29.6 4 13.5 4 29.6 4 13.5 29.6 29.6 ...
##  $ ORGANIC  : Factor w/ 3 levels "12","13.5","24.2": 2 1 3 2 3 1 3 2 1 1 ...
##  $ SITE_ID  : chr  "Ag2" "Ag1" "Ag3" "Ag2" ...
##  $ Kingdom  : chr  "Bacteria" "Bacteria" "Bacteria" "Bacteria" ...
##  $ Phylum   : chr  "Actinobacteria" "Actinobacteria" "Actinobacteria" "Proteobacteria" ...
##  $ Class    : chr  "Actinobacteria" "Actinobacteria" "Actinobacteria" "Alphaproteobacteria" ...
##  $ Order    : chr  "Streptomycetales" "Streptomycetales" "Streptomycetales" "Rhizobiales" ...
##  $ Family   : chr  "Streptomycetaceae" "Streptomycetaceae" "Streptomycetaceae" "Bradyrhizobiaceae" ...
##  $ Genus    : chr  "Streptomyces" "Streptomyces" "Streptomyces" "Bradyrhizobium" ...
# Relativos_S
# Conversión de valores Absolutos_S a Relativos_S
Transfomacion <- transform_sample_counts(Genero, function(x) x*100 / sum(x))
# Generamos el data.frame
V.Relativos_S <- psmelt(Transfomacion)
str(V.Relativos_S)
## 'data.frame':    3264 obs. of  16 variables:
##  $ OTU      : chr  "1883" "1883" "1883" "374" ...
##  $ Sample   : chr  "Ag2" "Ag3" "Ag1" "Ag2" ...
##  $ Abundance: num  11.1 9.93 9.83 7.4 6.96 ...
##  $ SAMPLE   : chr  "Ag2" "Ag3" "Ag1" "Ag2" ...
##  $ SITE     : chr  "3" "1" "2" "3" ...
##  $ pH       : num  6.7 6.1 6.2 6.7 6.1 6.2 6.1 6.2 6.2 6.1 ...
##  $ N_TOTAL  : num  0.22 0.29 0.36 0.22 0.29 0.36 0.29 0.36 0.36 0.29 ...
##  $ P_OLSEN  : num  13.5 4 29.6 13.5 4 29.6 4 29.6 29.6 4 ...
##  $ ORGANIC  : Factor w/ 3 levels "12","13.5","24.2": 2 3 1 2 3 1 3 1 1 3 ...
##  $ SITE_ID  : chr  "Ag2" "Ag3" "Ag1" "Ag2" ...
##  $ Kingdom  : chr  "Bacteria" "Bacteria" "Bacteria" "Bacteria" ...
##  $ Phylum   : chr  "Actinobacteria" "Actinobacteria" "Actinobacteria" "Proteobacteria" ...
##  $ Class    : chr  "Actinobacteria" "Actinobacteria" "Actinobacteria" "Alphaproteobacteria" ...
##  $ Order    : chr  "Streptomycetales" "Streptomycetales" "Streptomycetales" "Rhizobiales" ...
##  $ Family   : chr  "Streptomycetaceae" "Streptomycetaceae" "Streptomycetaceae" "Bradyrhizobiaceae" ...
##  $ Genus    : chr  "Streptomyces" "Streptomyces" "Streptomyces" "Bradyrhizobium" ...
# Renombramos géneros, phylum y familias que tengan menor al 1 %
V.Relativos_S$Genus[V.Relativos_S$Abundance < 1.0] <- "Genus < 1.0 % abund."
V.Relativos_S$Family[V.Relativos_S$Abundance < 1.0] <- "Family < 1.0 % abund."
V.Relativos_S$Phylum[V.Relativos_S$Abundance < 1.0] <- "Phyla < 1.0 % abund."

Graficamos la abundancia a diferente nivel taxonómico

Taxon <- c("Family", "Phylum")
# Ciclo
for (t in Taxon){
  # Prueba
  # t <- c("Family")
  # Absoluto
  plot.a <- ggplot(data = V.Absolutos_S, aes(x = V.Relativos_S[,t], y = Abundance, 
                                           fill = SITE_ID)) +
  geom_bar(aes(), stat = "identity", position = "stack") + 
  theme_classic() +
  labs(x="Phylum", y="Absolute abundance",
       title = "Absolute abundance of analyzed samples") +
  theme(legend.position = "right", 
        legend.text = element_text(size=12, face="bold"),
        legend.title = element_text(size = 12, face = "bold"),
        axis.text = element_text(color = "black", size = 12, face = "bold"),
        axis.title = element_text(color = "black", size = 12, face = "bold"),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        plot.title = element_text(hjust = 0.5, size=14, face="bold"))
  print(plot.a)
  pdf(paste0("Resultados/PlotAbunAbs",t,".pdf"), width=20, height=10)
  print(plot.a)
  dev.off()
  svg(paste0("Resultados/PlotAbunAbs",t,".svg"), width=20, height=10)
  print(plot.a)
  dev.off()
  # Relativo
  plot.r <- ggplot(data = V.Relativos_S, aes(x = V.Relativos_S[,t], y = Abundance, 
                                           fill = SITE_ID)) + 
    geom_bar(aes(), stat = "identity", position = "fill") + 
    theme_classic() + 
    labs(x="Phylum", y="Relative abundance",
         title = "Relative abundance of analyzed samples") +
    theme(legend.position = "right", 
      legend.text = element_text(size=12, face="bold"),
      legend.title = element_text(size = 12, face = "bold"),
      axis.text = element_text(color = "black", size = 12, face = "bold"),
      axis.title = element_text(color = "black", size = 12, face = "bold"),
      axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
      plot.title = element_text(hjust = 0.5, size=14, face="bold"))
  print(plot.r)
  pdf(paste0("Resultados/PlotAbunRel",t,".pdf"), width=20, height=10)
  print(plot.r)
  dev.off()
  svg(paste0("Resultados/PlotAbunRel",t,".svg"), width=20, height=10)
  print(plot.r)
  dev.off()
}
Gráficos de abundancia a diferentes niveles taxonómicos.

Gráficos de abundancia a diferentes niveles taxonómicos.

Gráficos de abundancia a diferentes niveles taxonómicos.

Gráficos de abundancia a diferentes niveles taxonómicos.

Gráficos de abundancia a diferentes niveles taxonómicos.

Gráficos de abundancia a diferentes niveles taxonómicos.

Gráficos de abundancia a diferentes niveles taxonómicos.

Gráficos de abundancia a diferentes niveles taxonómicos.

Análisis exploratorio de datos de microbiomas

Trabajamos a partir del objeto phyloseq

physeq <- merged_metagenomes2

Alfa diversidad

# Alfa deversidad
plot_richness(physeq, x = "SITE_ID", color = "SITE_ID")

Abundancias

# Abundancias
plot_bar(physeq, fill = "SITE_ID")

# Nos da el nombre de los taxones top 5
TopNGenus <- names(sort(taxa_sums(physeq), TRUE)[1:5])
# Generamos un objeto phyloseq con solo los 5 taxa que escogimos
Top5Genus <- prune_taxa(TopNGenus, physeq)
# Graficamos
plot_bar(Top5Genus, fill = "SITE_ID")

plot_bar(Top5Genus, fill = "SITE_ID", facet_grid = ~SITE_ID)

#Abundancias relativas
top_5_relativo <- transform_sample_counts(Top5Genus, function(x) x/sum(x)*100)
#Grafico de barras con abundancias relativas
plot_bar(top_5_relativo, fill="SITE_ID")

plot_bar(top_5_relativo, fill="SITE_ID", facet_grid = ~SITE_ID)

#Convertimos a dataframe
top_5_relativo_df <- psmelt(top_5_relativo)
#Damos formato como caracter y factor a los OTU
top_5_relativo_df$OTU <- as.character(top_5_relativo_df$OTU)
top_5_relativo_df$OTU <- as.factor(top_5_relativo_df$OTU)

#Convertimos a factor las muestras, no fue necesario determinar los niveles, ya los
#detecta como tal
top_5_relativo_df$Sample <- as.factor(top_5_relativo_df$Sample)

#definimos paleta de colores por niveles de OTU
OTU_colors<- brewer.pal(length(levels(top_5_relativo_df$OTU)),"Dark2")
#especificamos los aes de ggplot
global_rel_plot<-ggplot(top_5_relativo_df, aes(x=Sample, y=Abundance, fill=OTU)) +
  geom_bar(aes(), stat="identity", position="stack") +
  scale_fill_manual(values=OTU_colors) +
  labs(x= "Sample", y = "Relative abundance") +
  facet_grid(~SITE_ID, scales = "free_x") +
  ggtitle("Relative abundance by Groups") +
  theme(axis.title = element_text(size = 14), 
        axis.text.x = element_text(size = 5), 
        plot.title = element_text(size = 20))

global_rel_plot

Gráficos de mapa de calor

# Nos da el nombre de los taxones top 5
TopNGenus <- names(sort(taxa_sums(physeq), TRUE)[1:5])
# Generamos un objeto phyloseq con solo los 5 taxa que escogimos
Top5Genus <- prune_taxa(TopNGenus, physeq)
# Graficamos
plot_heatmap(Top5Genus)
## Warning in metaMDS(veganifyOTU(physeq), distance, ...): stress is (nearly) zero:
## you may have insufficient data
## Warning in postMDS(out$points, dis, plot = max(0, plot - 1), ...): skipping
## half-change scaling: too few points below threshold

plot_heatmap(top_5_relativo)
## Warning in metaMDS(veganifyOTU(physeq), distance, ...): stress is (nearly) zero:
## you may have insufficient data

## Warning in metaMDS(veganifyOTU(physeq), distance, ...): skipping half-change
## scaling: too few points below threshold

# Método NMDS y distancia de Bray-Curtis
p <- plot_heatmap(Top5Genus, "NMDS", "bray")
## Warning in metaMDS(veganifyOTU(physeq), distance, ...): stress is (nearly) zero:
## you may have insufficient data

## Warning in metaMDS(veganifyOTU(physeq), distance, ...): skipping half-change
## scaling: too few points below threshold
p

plot_heatmap(Top5Genus, "NMDS", "bray", low="#000033", high="#CCFF66")
## Warning in metaMDS(veganifyOTU(physeq), distance, ...): stress is (nearly) zero:
## you may have insufficient data

## Warning in metaMDS(veganifyOTU(physeq), distance, ...): skipping half-change
## scaling: too few points below threshold

plot_heatmap(Top5Genus, "NMDS", "bray", low="#000033", high="#FF3300")
## Warning in metaMDS(veganifyOTU(physeq), distance, ...): stress is (nearly) zero:
## you may have insufficient data

## Warning in metaMDS(veganifyOTU(physeq), distance, ...): skipping half-change
## scaling: too few points below threshold

plot_heatmap(Top5Genus, "NMDS", "bray", low="#000033", high="#66CCFF")
## Warning in metaMDS(veganifyOTU(physeq), distance, ...): stress is (nearly) zero:
## you may have insufficient data

## Warning in metaMDS(veganifyOTU(physeq), distance, ...): skipping half-change
## scaling: too few points below threshold

plot_heatmap(Top5Genus, "NMDS", "bray", low="#66CCFF", high="#000033", na.value="white")
## Warning in metaMDS(veganifyOTU(physeq), distance, ...): stress is (nearly) zero:
## you may have insufficient data

## Warning in metaMDS(veganifyOTU(physeq), distance, ...): skipping half-change
## scaling: too few points below threshold

# Método de PCoA y distancia de Bray-Curtis
plot_heatmap(Top5Genus, "PCoA", "bray")

Gráfico de redes

set.seed(123)
# Hacemos una gráfica a partir de el objeto phyloseq
ig <- make_network(physeq, max.dist = 0.55)
# Graficamos
plot_network(ig, physeq, color="SITE_ID", shape="SITE_ID")

# Graficamos sin necesidad de hacer el objeto anterior
plot_net(physeq, maxdist = .3, color = "SITE_ID", shape="SITE_ID")

Agrupamiento

abund_table <- t(physeq@otu_table@.Data)
abund_table_norm <- decostand(abund_table, "normalize")
bc_dist <- vegdist(abund_table_norm , method = "bray")

# Sencillo
cluster_single <- hclust (bc_dist, method = 'single')
plot(cluster_single)

# Completo
cluster_complete <- hclust (bc_dist, method = 'complete')
plot(cluster_complete)

# Promedio
cluster_average <- hclust (bc_dist, method = 'average')
plot(cluster_average)

# Mínima de Ward
cluster_ward <- hclust (bc_dist, method = 'ward.D2')
plot(cluster_ward)

# Todos juntos
# Correrlo en consola para visualizarlo mejor
par (mfrow = c(2,2))
plot(cluster_single)
plot(cluster_complete)
plot(cluster_average)
plot(cluster_ward)

par (mfrow = c(1,1))

Mismos gráficos usando el paquete factorextra Mapa de calor

fviz_dist(dist.obj = bc_dist, lab_size = 8)

set.seed(12345)

Dendogramas

# Sencilla
hc_single <- hclust(d=bc_dist, method = "single")
fviz_dend(x = hc_single, k=3,
cex = 0.7,
main = "",
xlab = "Samples",
ylab = "Distance",
sub = "",
horiz = TRUE)+
geom_hline(yintercept = 0.5, linetype = "dashed")
## Registered S3 method overwritten by 'dendextend':
##   method     from 
##   rev.hclust vegan
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <]8;;https://github.com/kassambara/factoextra/issueshttps://github.com/kassambara/factoextra/issues]8;;>.

# Completa
hc_complete <- hclust(d=bc_dist, method = "complete")
fviz_dend(x = hc_complete, k=3,
cex = 0.7,
main = "",
xlab = "Samples",
ylab = "Distance",
sub = "",horiz = TRUE)+
geom_hline(yintercept = 0.55, linetype = "dashed")

# Ward
hc_ward <- hclust(d=bc_dist, method = "ward.D2")
fviz_dend(x = hc_ward, k=3,
cex = 0.7,
main = "",
xlab = "Samples",
ylab = "Distance",
sub = "",
horiz = TRUE)+
geom_hline(yintercept = 0.55, linetype = "dashed")

# Promedio
hc_average <- hclust(d=bc_dist, method = "average")
fviz_dend(x = hc_average, k=3,
cex = 0.7,
main = "",
xlab = "Samples",
ylab = "Distance",
sub = "",
horiz = TRUE) +
geom_hline(yintercept = 0.55, linetype = "dashed")

Ordenación

Lo vamos aplicar para todos los OTUs y para los Top 10

Top10 <- names(sort(taxa_sums(physeq), TRUE)[1:10])
# Generamos un objeto phyloseq con solo los 5 taxa que escogimos
Top10Genus <- prune_taxa(Top10, physeq)
Abund10 <- Top10Genus@otu_table@.Data

Análisis de Componentes Principales

Todos

# Tabla de metadatos
meta_table <- MetaShotgun
# Estandarizamos la abundancia de los datos
stand_abund_table <- decostand(abund_table, method = "total")
PCA <- rda(stand_abund_table)
PCA
## Call: rda(X = stand_abund_table)
## 
##                 Inertia Rank
## Total         0.0001378     
## Unconstrained 0.0001378    2
## Inertia is variance 
## 
## Eigenvalues for unconstrained axes:
##        PC1        PC2 
## 0.00012967 0.00000818
# Con el siguiente código podemos pedir los valores de los eigenvalores.
eig <- PCA$CA$eig
eig
##            PC1            PC2 
## 0.000129668682 0.000008175732
# Variación total
sum(apply(stand_abund_table, 2, var))
## [1] 0.0001378444
# Otra manera de visualizar
head(summary(PCA))
## 
## Call:
## rda(X = stand_abund_table) 
## 
## Partitioning of variance:
##                 Inertia Proportion
## Total         0.0001378          1
## Unconstrained 0.0001378          1
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                             PC1         PC2
## Eigenvalue            0.0001297 0.000008176
## Proportion Explained  0.9406887 0.059311306
## Cumulative Proportion 0.9406887 1.000000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  0.128856 
## 
## 
## Species scores
## 
##                PC1       PC2
## 374      0.0325550 0.0049411
## 1076    -0.0069872 0.0021648
## 1842539 -0.0012303 0.0005328
## 1882747 -0.0013264 0.0002074
## 912     -0.0013639 0.0004397
## 1333996 -0.0004035 0.0006266
## ....                        
## 
## 
## Site scores (weighted sums of species scores)
## 
##          PC1       PC2
## Ag1 -0.05976 -0.086593
## Ag2  0.10487 -0.008455
## Ag3 -0.04511  0.095048
# variables (species) tienen la mayor correlación absoluta al primer y segundo eje:
head(sort(abs(PCA$CA$v[,1]),decreasing = TRUE))
##      1883      1386    191495       374    196162      1873 
## 0.6208206 0.3542973 0.3153027 0.2604890 0.2599363 0.1916209
head(sort(abs(PCA$CA$v[,2]),decreasing = TRUE))
##     36814      1763      1873      1386       286    331696 
## 0.4962678 0.4471708 0.3458962 0.2481418 0.2158779 0.1865143
# Diagramas utilizando la función biplot(). 
# La opción de visualización “species” es la etiqueta del paquete vegan para OTUs/taxa. 
# Por defecto es “sites” (etiqueta para muestras).
biplot(PCA, display = 'species', type = c("text"))
## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

#Mismo biplot con etiquetas
biplot(PCA, display = c('sites','species'), type=c("text","points"))
## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

## Warning in arrows(0, 0, g$species[, 1] * arrlen, g$species[, 2] * arrlen, :
## zero-length arrow is of indeterminate angle and so skipped

ordiplot(PCA, display = "sites", type = "text")

Top10

# Tabla de metadatos
meta_table <- MetaShotgun
# Estandarizamos la abundancia de los datos
stand_abund_table <- decostand(Abund10, method = "total")
PCA <- rda(stand_abund_table)
PCA
## Call: rda(X = stand_abund_table)
## 
##                Inertia Rank
## Total         0.002636     
## Unconstrained 0.002636    2
## Inertia is variance 
## 
## Eigenvalues for unconstrained axes:
##       PC1       PC2 
## 0.0024283 0.0002073
# Con el siguiente código podemos pedir los valores de los eigenvalores.
eig <- PCA$CA$eig
eig
##          PC1          PC2 
## 0.0024283334 0.0002072756
# Variación total
sum(apply(stand_abund_table, 2, var))
## [1] 0.002635609
# Otra manera de visualizar
head(summary(PCA))
## 
## Call:
## rda(X = stand_abund_table) 
## 
## Partitioning of variance:
##                Inertia Proportion
## Total         0.002636          1
## Unconstrained 0.002636          1
## 
## Eigenvalues, and their contribution to the variance 
## 
## Importance of components:
##                            PC1       PC2
## Eigenvalue            0.002428 0.0002073
## Proportion Explained  0.921356 0.0786443
## Cumulative Proportion 0.921356 1.0000000
## 
## Scaling 2 for species and site scores
## * Species are scaled proportional to eigenvalues
## * Sites are unscaled: weighted dispersion equal on all dimensions
## * General scaling constant of scores:  0.3924469 
## 
## 
## Species scores
## 
##         PC1       PC2
## Ag1 -0.1656 -0.075724
## Ag2  0.3073 -0.004037
## Ag3 -0.1417  0.079762
## 
## 
## Site scores (weighted sums of species scores)
## 
##              PC1       PC2
## 374     0.001874 -0.038646
## 407    -0.095652 -0.010277
## 32008  -0.078443 -0.142850
## 286    -0.112012 -0.180988
## 1883    0.036769 -0.064459
## 1763   -0.002113  0.176748
## 36814  -0.020847  0.234887
## 1873   -0.168049  0.078471
## 196162  0.219928 -0.002675
## 191495  0.218545 -0.050211
# variables (species) tienen la mayor correlación absoluta al primer y segundo eje:
head(sort(abs(PCA$CA$v[,1]),decreasing = TRUE))
##       Ag2       Ag1       Ag3 
## 0.8156721 0.4396057 0.3760663
head(sort(abs(PCA$CA$v[,2]),decreasing = TRUE))
##        Ag3        Ag1        Ag2 
## 0.72473498 0.68805047 0.03668451
# Diagramas utilizando la función biplot(). 
# La opción de visualización “species” es la etiqueta del paquete vegan para OTUs/taxa. 
# Por defecto es “sites” (etiqueta para muestras).
biplot(PCA, display = 'species', type = c("text"))

#Mismo biplot con etiquetas
biplot(PCA, display = c('sites','species'), type=c("text","text"))

ordiplot(PCA, display = "sites", type = "text")

Análisis de coordenadas principales (PCoA)

Todos

bc_dist <-vegdist(abund_table, "bray")
PCoA <- cmdscale(bc_dist, eig = TRUE, k = 2)
PCoA
## $points
##            [,1]          [,2]
## Ag1  0.03048085  0.0118099989
## Ag2 -0.05615181  0.0006942406
## Ag3  0.02567097 -0.0125042395
## 
## $eig
## [1]  4.741107e-03  2.963141e-04 -3.023333e-19
## 
## $x
## NULL
## 
## $ac
## [1] 0
## 
## $GOF
## [1] 1 1
# Variación explicada primer eje
explainedvar1 <- round(PCoA$eig[1] / sum(PCoA$eig), 2) * 100
explainedvar1
## [1] 94
# Variación explicada segundo eje
explainedvar2 <- round(PCoA$eig[2] / sum(PCoA$eig), 2) * 100
explainedvar2
## [1] 6
# Suma
sum_eig <- sum(explainedvar1, explainedvar2)
sum_eig
## [1] 100

Criterios de evaluación

# Correr todo junto
# Definir los parametros del plot
par(mar = c(5, 5, 1, 2) + 0.1)
# Fraficar eigenvalores
plot(PCoA$eig, xlab = "PCoA", ylab = "Eigenvalue",
las = 1, cex.lab = 1.5, pch = 16)
# Añadir expextación del criterio de Kaiser-Guttman y el modelo Broken Stick
abline(h = mean(PCoA$eig), lty = 2, lwd = 2, col = "blue")
b_stick <- bstick(8, sum(PCoA$eig))
lines(1:8, b_stick, type = "l", lty = 4, lwd = 2, col = "red")
# Añadir legendas
legend("topright", legend = c("Avg Eigenvalue", "Broken-Stick"),
lty = c(2, 4), bty = "n", col = c("blue", "red"))

Gráfico

# Correr todo junto
# Definir parámetros
par(mar = c(5, 5, 1, 2) + 0.1)
# Initiate Plot
plot(PCoA$points[ ,1], PCoA$points[ ,2], ylim = c(-0.5, 0.5),
xlab = paste("PCoA 1 (", explainedvar1, " %)", sep = ""),
ylab = paste("PCoA 2 (", explainedvar2, " %)", sep = ""),
pch = 5, cex = 1.0, type = "n", cex.lab = 1.0, cex.axis = 1.2, axes = FALSE)
# Añadir ejes
axis(side = 1, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
axis(side = 2, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
abline(h = 0, v = 0, lty = 3)
box(lwd = 2)
# Añadir puntos y texto
points(PCoA$points[ ,1], PCoA$points[ ,2],
pch = 19, cex = 3, bg = "blue", col = "blue")
text(PCoA$points[ ,1], PCoA$points[ ,2],
labels = row.names(PCoA$points))

# Correr todo junto
# Definir parámetros
par(mar = c(5, 5, 1, 2) + 0.1)
# Inicializar el plot
plot(PCoA$points[ ,1], PCoA$points[ ,2], ylim = c(-0.5, 0.5),
xlab = paste("PCoA 1 (", explainedvar1, " %)", sep = ""),
ylab = paste("PCoA 2 (", explainedvar2, " %)", sep = ""),
pch = 5, cex = 1.0, type = "n", cex.lab = 1.0, cex.axis = 1.2, axes = FALSE)
# Añadir ejes
axis(side = 1, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
axis(side = 2, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
abline(h = 0, v = 0, lty = 3)
box(lwd = 2)
# Añadir puntos y etiquetas
points(PCoA$points[ ,1], PCoA$points[ ,2],
pch = 19, cex = 3, bg = "blue", col = "blue")
text(PCoA$points[ ,1], PCoA$points[ ,2],
labels = row.names(PCoA$points))
CalakmulREL <- abund_table
for(i in 1:nrow(abund_table)){
CalakmulREL[i, ] = abund_table[i, ] / sum(abund_table[i, ])
}
# Calcular y añadir el score de las especies
PCoA <- add.spec.scores(PCoA,CalakmulREL,method = "pcoa.scores",Rscale=TRUE,scaling=1, multi=1)
text(PCoA$cproj[ ,1], PCoA $cproj[ ,2],
labels = row.names(PCoA$cproj),cex=0.5, col = "blue")

genus_corr <- add.spec.scores(PCoA, CalakmulREL, method = "cor.scores")$cproj
corrcut <- 0.8 #definir el cutoff  
import_genus <- genus_corr[abs(genus_corr[, 1]) >= corrcut | abs(genus_corr[, 2]) >= corrcut, ]
# Los 12 géneros importantes con correlación mayor o igual a 0,8 
# a lo largo de los ejes PCoA se imprimen a continuación:
import_genus[complete.cases(import_genus),] %>%
kbl(booktabs = TRUE, align = "c",
caption = "Géneros más importantes con correlación mayor o igual a 0.8" ) %>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 7)
Géneros más importantes con correlación mayor o igual a 0.8
Dim1 Dim2
374 -0.9835451 -0.1806627
1076 0.9639275 -0.2661649
1842539 0.9295353 -0.3687332
1882747 0.9923099 -0.1237786
912 0.9608191 -0.2771763
1333996 0.5672543 -0.8235427
379 -0.9998699 -0.0161320
358 0.8920150 -0.4520059
399 -0.8576021 0.5143138
380 0.9803048 -0.1974905
879274 0.8963447 -0.4433579
1273132 0.0484755 -0.9988244
407 0.9115661 -0.4111536
1882682 -0.9995394 -0.0303465
223967 0.9388086 -0.3444392
674703 -0.8948257 -0.4464157
53399 0.0867456 -0.9962305
1736675 0.5238160 -0.8518315
1069 0.9550470 -0.2964543
531813 0.9996074 -0.0280188
1608628 0.9883760 -0.1520293
454601 0.9301483 0.3671840
68287 -0.9983022 -0.0582463
83263 -0.3656077 -0.9307690
1867719 0.4969819 -0.8677609
1620421 0.9982966 -0.0583436
266779 0.9986729 -0.0515015
331696 0.8798245 -0.4752986
921 0.9943553 0.1061013
280 0.9120910 -0.4099878
7 0.9039212 -0.4276989
2170729 -0.0663632 -0.9977955
187303 0.9251290 -0.3796530
426 0.9799034 -0.1994727
1702325 0.9454020 -0.3259066
293089 0.9846307 0.1746491
256618 0.9994781 0.0323030
1922226 0.9976734 0.0681753
1458461 0.8120099 -0.5836436
2304600 0.9988035 -0.0489044
1798205 0.9974045 -0.0720018
199596 0.3739416 -0.9274523
533 0.8740029 -0.4859208
529 0.5233525 -0.8521163
234 0.0706983 -0.9974977
1686310 0.8641238 -0.5032794
160791 -0.9868568 -0.1615975
165697 -0.8834057 -0.4686090
165695 -0.9164974 -0.4000406
158500 -0.9996369 -0.0269447
1850238 -0.9137016 -0.4063858
150203 -0.8549144 0.5187691
2077182 -0.9706375 -0.2405471
1634516 -0.9688452 -0.2476671
692370 -0.9931837 -0.1165595
2003315 -0.8372654 -0.5467967
450378 -0.9045881 0.4262867
191 0.9706196 0.2406191
55518 0.9827923 0.1847143
34018 0.9759914 0.2178089
171437 0.9212596 0.3889483
28077 0.9990421 0.0437590
2220096 0.9988427 0.0480972
1612173 0.9934958 0.1138683
1084 0.9884813 -0.1513431
1288970 0.9752221 -0.2212282
220697 0.9240809 0.3821970
1549855 0.8037673 0.5949439
1263979 0.8542575 -0.5198501
257708 0.8272769 0.5617944
522 0.9504749 0.3108014
1434011 0.9437466 0.3306697
1969806 0.9290266 0.3700128
435 0.9417183 0.3364026
364410 0.9037383 0.4280854
320497 0.1745915 0.9846409
153496 0.9988196 0.0485745
91915 0.9353160 0.3538134
266 0.9800025 0.1989852
1063 0.9782096 -0.2076198
1208324 0.9566259 0.2913194
1917485 0.9487622 0.3159909
311180 0.9137252 -0.4063328
2021862 0.8341500 -0.5515376
60890 0.9914607 0.1304056
121719 0.9514046 -0.3079436
1229727 0.5204761 -0.8538762
1881061 0.9960067 -0.0892782
89184 0.9988127 -0.0487147
2009329 0.9663096 -0.2573824
1250539 0.9989293 -0.0462633
42444 0.9159171 -0.4013675
1335048 0.9631806 0.2688552
1915078 0.8799196 0.4751226
2109625 0.8687051 -0.4953296
1267768 0.9756921 -0.2191457
1609966 0.5924697 0.8055927
92945 0.8695219 -0.4938945
159346 0.3232299 0.9463205
2169400 -0.4507368 0.8926569
133924 0.9182266 -0.3960553
290400 0.9000987 -0.4356861
299262 0.9981937 0.0600782
74030 0.9076401 0.4197493
2434 0.9321858 -0.3619803
53946 -0.8976650 0.4406785
245188 0.9999002 0.0141266
379347 -0.0981275 0.9951739
441209 0.1589096 -0.9872931
1920883 0.2479234 -0.9687796
1873716 0.8821695 -0.4709321
74318 0.9076378 -0.4197542
2724 0.9995238 0.0308586
69666 -0.0741974 -0.9972436
588932 0.3939680 -0.9191241
284016 0.1078947 -0.9941623
78587 0.4528011 -0.8916116
1868589 0.9455983 0.3253365
991904 0.9551367 -0.2961653
349221 0.4844556 -0.8748158
767892 0.8600871 0.5101472
953 0.9787591 -0.2050136
779 -0.8036887 0.5950499
33994 0.5443255 0.8388741
142058 -0.9948845 0.1010186
784 -0.9823412 0.1870982
2115978 0.3072179 -0.9516392
1802984 0.9611989 -0.2758563
1528098 0.9511609 -0.3086956
208216 0.9838229 -0.1791436
1124597 0.9983725 0.0570291
744985 -0.3294819 0.9441619
91604 0.5903187 0.8071703
1509244 0.5825715 -0.8127795
164546 0.9779951 0.2086277
1822464 0.9656840 0.2597199
48736 0.9856384 0.1688695
93220 0.9538928 0.3001476
576610 0.9984642 -0.0554008
29443 0.9997936 -0.0203150
2268024 0.3655372 0.9307967
34073 0.9992231 0.0394109
80867 0.9999167 0.0129075
94132 0.9993583 -0.0358190
795665 0.9702043 0.2422884
1842727 0.9991048 -0.0423039
285 0.9999211 0.0125586
80865 0.8495420 -0.5275211
296591 0.9999927 0.0038341
2109913 0.9962544 -0.0864707
179636 0.9202464 -0.3913395
364317 0.9259870 0.3775554
2109914 0.9820725 0.1885035
2109915 0.9777052 -0.2099821
1678128 0.9943941 0.1057369
1844971 0.8232774 0.5676393
1436290 0.9900782 0.1405177
1546149 0.8754699 0.4832726
517 0.9800336 0.1988322
85698 0.9742677 0.2253941
75697 0.9988167 -0.0486325
1851544 0.9689941 0.2470837
302406 0.9854930 0.1697161
1007105 0.9833795 0.1815621
90245 0.8491953 -0.5280789
1472345 0.1872288 0.9823163
29575 0.8278955 -0.5608824
28068 0.9687759 -0.2479381
946333 0.9809114 -0.1944550
34029 0.9998204 -0.0189542
1050370 0.9964469 0.0842234
105560 0.9866938 0.1625894
1658665 0.9376489 0.3475840
1296669 0.8143521 0.5803712
76731 0.9675133 -0.2528200
864051 0.9947296 0.1025328
1141883 0.9357519 0.3526589
92645 -0.3520820 0.9359692
55508 0.8433815 -0.5373152
158899 0.9476497 0.3193118
847 0.9996849 -0.0251020
2259134 0.9873975 0.1582595
748247 0.9998373 0.0180356
2005884 0.9759244 0.2181089
551760 0.9697469 0.2441125
146939 0.9983865 -0.0567834
259537 0.9996583 -0.0261401
535 0.9938822 -0.1104449
1906741 0.8811667 -0.4728057
748280 0.9920682 0.1257011
57480 0.9999939 0.0034860
1192162 0.9325883 -0.3609420
168471 0.9201738 -0.3915101
1938604 0.4793148 0.8776430
332411 0.9623251 -0.2719013
2290923 0.9053492 0.4246680
482 0.9633857 0.2681195
63 0.9269176 -0.3752649
539 0.9853479 -0.1705568
72 0.9061647 -0.4229251
748811 0.9667525 -0.2557139
1842540 0.9807147 -0.1954447
63745 0.9388060 0.3444462
649841 0.9799198 0.1993921
1985873 0.9997699 -0.0214527
1188319 0.9992178 0.0395450
370405 0.9919357 0.1267418
36861 0.9634981 0.2677152
44574 0.9966370 0.0819433
1288494 0.9961584 0.0875695
81682 0.9449092 -0.3273326
1662285 0.8514902 0.5243705
359408 0.9897576 -0.1427582
1904640 0.9997951 0.0202439
327160 0.9548926 0.2969513
33056 0.9710479 0.2388849
189385 0.9970087 -0.0772896
1160784 -0.9078740 0.4192431
669502 0.4740952 -0.8804736
286 0.8640389 0.5034250
353 0.9980421 -0.0625452
1697053 -0.9914017 -0.1308536
1879049 0.4842060 -0.8749540
497 0.9974717 0.0710653
40324 0.0810075 0.9967135
69 0.5873723 -0.8093169
338 0.9985045 0.0546688
314722 0.8891440 0.4576276
445710 0.9325029 0.3611624
323415 -0.9551513 0.2961181
666685 0.3554083 -0.9347111
242606 0.3455498 0.9384004
2021234 0.8741273 0.4856968
550 -0.9968684 0.0790787
570 -0.8825882 0.4701468
413496 -0.8244764 0.5658963
544 -0.9678216 0.2516371
203907 0.9940635 0.1088014
1410383 0.4861943 0.8738507
1778262 -0.0663593 0.9977958
1778264 0.2491447 -0.9684663
1778263 -0.9486368 -0.3163671
472834 -0.9078740 0.4192431
208223 -0.4468888 0.8945895
28901 -0.9991548 -0.0411055
61646 -0.9892995 0.1458991
158822 -0.1716400 0.9851597
929813 -0.1364650 0.9906449
83654 -0.9210432 -0.3894605
563 0.5216857 0.8531378
2479367 0.3271624 -0.9449681
2172103 0.9739395 0.2268078
622 -0.2830481 0.9591057
1048758 -0.9597533 0.2808444
168169 0.9099969 0.4146151
613 0.9979288 0.0643288
629 -0.9085305 0.4178185
1805933 -0.9540848 0.2995366
1639108 0.9695117 -0.2450452
1878942 -0.4873568 -0.8732029
554 -0.1274744 0.9918419
1239307 -0.3297449 0.9440701
598467 0.0999721 -0.9949902
1082704 -0.9825653 0.1859180
470934 -0.9988709 0.0475075
182337 -0.9344836 0.3560062
665914 -0.8348387 0.5504947
82987 -0.9888229 0.1490949
51229 0.2471352 0.9689810
82983 0.9556418 -0.2945315
40576 0.1928302 0.9812321
587 -0.5125625 -0.8586499
582 -0.9244155 0.3813869
230089 -0.3985315 -0.9171546
584 -0.9939415 0.1099107
1972431 0.8157587 0.5783923
158841 0.8972662 0.4414900
108010 0.9847271 0.1741047
160660 0.9883339 0.1523025
351052 0.8755321 0.4831600
1335757 0.9805336 0.1963514
1442136 0.9980133 -0.0630041
1166950 0.9826809 0.1853057
80679 0.8892097 0.4574998
37487 0.8160903 0.5779244
1049 0.9898544 0.1420857
73141 0.9448672 0.3274539
133539 0.8641612 0.5032152
1630141 0.5455714 0.8380644
1579979 0.9382242 0.3460280
1860122 0.9895352 -0.1442918
1548547 0.9758654 0.2183730
475662 0.8727041 0.4882494
376489 0.4438329 0.8961096
1771309 0.9732784 -0.2296282
698828 0.8474523 0.5308716
28258 0.8077306 0.5895518
33074 -0.8278449 0.5609570
1495769 -0.5529811 0.8331938
2014542 0.8610438 0.5085307
1917421 0.9047439 0.4259560
178399 0.9875946 0.1570252
1249553 0.9684922 -0.2490440
187493 -0.5237264 0.8518866
2161747 0.9619789 -0.2731239
1336806 0.9992560 -0.0385662
1445505 0.3862156 0.9224085
158327 0.8727484 0.4881703
1561924 0.9664096 -0.2570068
1027273 0.4594553 0.8882009
642 0.8450593 0.5346726
347534 -0.0768462 -0.9970430
511062 0.8033403 0.5955203
1416627 -0.8379880 0.5456887
465721 0.9699204 0.2434221
2303331 0.9824620 0.1864628
1675686 0.9805973 0.1960328
1620215 0.9982788 -0.0586461
1281578 0.9952876 0.0969672
39775 0.9976549 0.0684449
416 0.9973673 -0.0725155
414 0.9303535 0.3666639
1432792 0.9986538 -0.0518712
1704499 0.9882672 0.1527348
1420917 0.9449006 0.3273574
28108 0.9658373 0.2591493
2183582 0.5242436 0.8515683
1526571 0.9384723 0.3453545
300231 0.9334532 -0.3586992
326544 -0.9457443 0.3249118
680279 -0.9264833 -0.3763359
359303 0.9632179 0.2687215
53246 0.5166018 -0.8562258
44012 0.9903573 -0.1385371
2100422 0.9913529 0.1312226
1763536 0.1751189 -0.9845473
1248727 0.9388135 0.3444259
585455 0.9368641 0.3496936
1543721 0.9966152 0.0822081
1076588 0.9985717 0.0534280
1249552 0.9273032 -0.3743111
2291597 -0.9384068 -0.3455326
1427364 -0.4890659 -0.8722468
186490 0.8688363 -0.4950994
1769779 0.9559117 0.2936543
1987723 0.9930523 0.1176740
447471 -0.9995202 -0.0309732
2426 0.9455627 0.3254399
1737490 0.5814416 0.8135882
930806 0.8674989 -0.4974392
393662 0.9152976 -0.4027782
1470434 -0.2123446 -0.9771948
1620392 -0.4458136 -0.8951258
716816 -0.3873631 -0.9219272
672 0.9966610 0.0816505
668 0.9997414 -0.0227397
673 0.1845721 0.9828190
1927128 -0.1089770 0.9940443
1755811 0.3818304 -0.9242324
1810504 0.8600261 0.5102500
2183911 0.9917927 -0.1278561
446 -0.4593416 -0.8882597
451 -0.1147652 -0.9933926
463 -0.9125532 -0.4089580
777 0.8063812 0.5913962
727 0.3775727 -0.9259800
67854 -0.5349491 0.8448843
75985 0.4338758 0.9009727
747 -0.2335891 -0.9723354
157673 0.8195287 0.5730381
731 0.9974510 -0.0713554
738 -0.9998697 0.0161431
47735 0.4521939 -0.8919197
754476 0.9962920 0.0860359
34067 0.5582005 -0.8297061
39765 0.0864185 -0.9962589
262 -0.9535931 -0.3010983
594679 0.9715788 0.2367162
288004 0.4852773 -0.8743603
40754 0.9802414 -0.1978049
1196095 -0.9998096 0.0195115
1267021 0.8679429 0.4966640
870 0.8657035 -0.5005572
404589 0.9999000 -0.0141409
43 0.8331507 0.5530460
41 0.8761220 0.4820894
83453 0.8917949 0.4524398
32 0.8667575 0.4987299
184914 0.8996916 0.4365260
1391653 0.8562075 0.5166320
56 0.9942470 0.1071116
52 0.9216759 0.3879607
1882918 0.9849700 0.1727253
927083 0.9742840 0.2253233
80816 0.9641180 0.2654741
443143 0.9904408 0.1379386
483547 0.9022484 0.4312167
1823759 0.8563186 0.5164479
29543 0.9940436 0.1089827
881 0.9876905 0.1564206
1716143 0.8075027 0.5898638
29546 0.0191292 -0.9998170
45663 0.9905909 0.1368562
897 0.9986291 -0.0523444
259354 0.9894677 -0.1447539
2296 0.8823439 0.4706053
2293 0.3641637 -0.9313350
28223 0.9999681 -0.0079878
427923 0.9981301 0.0611256
1986146 0.3477060 0.9376036
84980 0.0085180 -0.9999637
119484 0.8476997 0.5304765
2358 0.9904994 0.1375173
60893 0.9940087 0.1093008
316277 0.9219588 -0.3872880
453230 0.9422290 0.3349694
1548548 0.9854470 -0.1699829
1621989 0.9386684 0.3448212
33002 0.9861858 0.1656428
84405 0.5948440 -0.8038411
673862 -0.8836440 -0.4681594
28198 0.8798679 0.4752183
366522 0.5388493 -0.8424022
210 0.9707893 0.2399337
1176482 0.9893160 0.1457873
148813 -0.5604281 0.8282031
387093 -0.8808971 0.4733078
387092 0.1390396 0.9902868
244787 0.9871942 -0.1595231
959 -0.0334129 0.9994416
2109558 -0.0129370 0.9999163
1921087 0.9649601 0.2623965
1883 -0.9998396 -0.0179116
2126346 0.0880403 -0.9961169
1763 -0.5416584 -0.8405987
36814 -0.3125967 -0.9498859
36809 -0.9927541 0.1201634
875328 0.3395642 -0.9405829
639313 -0.9756662 -0.2192610
1827 -0.9993684 0.0355367
37326 -0.2204426 0.9754000
84595 0.3091210 -0.9510227
1528099 0.9981351 0.0610442
33882 -0.9733996 -0.2291140
589382 -0.9687273 -0.2481279
1575 -0.8998519 -0.4361956
28447 0.9688878 0.2475005
399736 -0.9215517 0.3882556
69373 -0.5614267 -0.8275264
1795630 0.3179108 -0.9481206
1978566 -0.9277831 0.3731200
412690 0.0400547 -0.9991975
2419774 -0.4897283 -0.8718751
2419771 -0.9994046 -0.0345019
150123 -0.4756757 -0.8796207
33888 -0.8328076 -0.5535625
279828 0.4470229 -0.8945225
1619308 -0.9491993 0.3146754
2079791 0.8157703 -0.5783760
1987356 -0.9970123 0.0772425
535712 0.2686331 -0.9632426
1855377 0.9969458 -0.0780961
1159327 0.9958073 0.0914754
1357915 0.9148126 -0.4038786
446860 0.9027098 0.4302499
121292 0.5376246 -0.8431843
37927 0.8787472 -0.4772874
37929 0.9970826 0.0763300
2058657 -0.9494746 -0.3138438
43675 0.2434047 -0.9699248
43663 0.9790037 0.2038427
2170745 0.3795983 0.9251514
1618207 -0.9874253 0.1580860
1646 0.1649868 0.9862958
53358 0.9961769 -0.0873587
2283195 0.9953364 0.0964644
1658671 0.3392316 -0.9407029
186189 -0.3291975 -0.9442611
1708 -0.8815753 -0.4720434
545619 0.5322578 -0.8465823
1903186 0.9999053 -0.0137617
1630135 0.9734695 0.2288170
472569 0.9650343 -0.2621238
571913 -0.9897269 0.1429710
1274 0.9047572 -0.4259276
1276 0.9558675 0.2937981
1863 -0.4631708 0.8862691
43674 0.5477096 -0.8366685
2039 0.9972159 -0.0745687
1873 0.9231768 -0.3843755
1865 0.1849835 0.9827416
673534 0.9737269 0.2277190
168697 0.9626805 0.2706404
1003110 0.9666411 -0.2561346
33910 -0.5838677 0.8118488
240495 -0.9935899 0.1130447
2072503 -0.9517815 0.3067766
40566 -0.9727667 0.2317866
860235 -0.8849047 0.4657722
103731 -0.9621625 0.2724762
43357 -0.3356927 0.9419715
1836 -0.9989485 0.0458462
211114 -0.9598572 0.2804890
1653480 -0.3080125 0.9513823
196162 -0.9972033 -0.0747364
2079793 -0.9823230 -0.1871936
2045 -0.9935856 -0.1130826
117157 0.9719907 -0.2350191
642780 -0.9981142 -0.0613842
182640 -0.9907254 0.1358792
546871 -0.9989776 0.0452070
75385 0.5218918 -0.8530117
1909732 0.9931842 -0.1165557
630515 -0.9994133 -0.0342485
33010 0.8521848 0.5232409
1744 0.8168815 -0.5768056
1750 0.1540603 -0.9880614
1871034 0.9659141 0.2588628
2202249 0.9618120 -0.2737109
1909395 0.3690083 0.9294261
1411117 0.0207671 0.9997843
2020 0.4915804 0.8708322
2014 0.8809777 0.4731578
298654 0.9958381 -0.0911399
1907575 -0.9669613 -0.2549230
1861 0.9893206 -0.1457557
477641 0.9065652 -0.4220657
138336 0.9995656 -0.0294714
419479 -0.2133579 -0.9769741
111015 0.9544000 0.2985309
59505 0.9868178 -0.1618349
28264 0.3415460 0.9398651
2051 0.3393901 -0.9406457
53461 0.1519176 0.9883932
1678 0.9950613 0.0992620
78258 0.9696634 -0.2444441
2702 0.9256662 -0.3783412
78259 -0.1894001 -0.9819000
2006 0.9705884 -0.2407451
283811 0.4033146 0.9150614
131568 0.4725685 -0.8812940
414996 -0.5175073 0.8556788
1884913 0.4178620 0.9085105
2039639 0.5502179 0.8350211
191495 -0.9998097 -0.0195101
1841863 0.9351239 -0.3543209
84112 0.9960247 0.0890778
446660 -0.9042642 0.4269734
84110 0.9921057 -0.1254044
79604 -0.9898497 0.1421183
1944646 -0.0230226 -0.9997349
1805478 0.9843119 -0.1764372
1871022 -0.4064108 -0.9136905
1382 -0.1436886 0.9896229
33871 -0.9449382 -0.3272488
74426 0.2314749 -0.9728409
49319 -0.9947905 -0.1019402
53635 0.8715400 -0.4903245
1386 -0.9897854 0.1425653
1221500 0.9999290 0.0119158
129337 0.9137092 0.4063687
1421 -0.2239581 0.9745988
163877 0.9095561 0.4155812
182710 0.8771414 0.4802322
1295642 0.8888178 0.4582608
586416 0.4354077 0.9002334
1230341 0.9791471 0.2031528
1570 -0.8888319 -0.4582334
33936 0.9991152 0.0420564
1472767 0.9973853 -0.0722678
1449 -0.3371328 0.9414571
2213194 -0.8740578 -0.4858220
61624 0.9814381 0.1917791
377615 0.9701233 0.2426126
1393 0.9149874 0.4034823
1450761 0.9126745 0.4086872
1903704 0.0591461 0.9982493
33943 0.9729464 0.2310310
1038856 0.8368894 0.5473721
1571 0.9336218 0.3582600
1508404 0.8390240 0.5440944
241244 0.9897344 0.1429190
1750719 -0.9788218 0.2047140
1280 0.9015666 0.4326404
1461582 0.4910534 0.8711295
407035 -0.9370821 -0.3491091
1849491 0.9897288 -0.1429577
360911 0.9925974 0.1214510
29391 0.3397568 0.9405133
1471761 0.9632862 -0.2684765
1639 0.8971603 0.4417050
2756 0.9347673 0.3552607
1584 0.8456507 0.5337367
1358 0.4881107 0.8727817
1352 0.9769710 0.2133719
33970 0.3584591 0.9335454
633807 -0.4726054 0.8812742
1243 0.2174864 0.9760633
137591 0.9897288 -0.1429577
2036206 0.0612856 0.9981203
1903686 0.8249070 0.5652685
1911586 0.9698597 0.2436640
1485 0.9900972 0.1403832
2086584 0.8175992 0.5757878
1424294 0.9173160 0.3981600
471827 0.9903224 0.1387857
1564 0.9621472 0.2725302
339862 0.9039422 0.4276547
49338 0.2745889 -0.9615617
102134 0.9940717 -0.1087265
56112 0.9913136 0.1315194
51197 0.9729735 0.2309169
58138 0.8003018 0.5995974
51515 0.8158837 0.5782160
863643 0.8778999 -0.4788443
73919 0.9460400 0.3240500
178899 0.8696400 0.4936865
86332 0.0152207 0.9998842
2734 0.9987992 0.0489918
292800 0.9977862 0.0665037
853 0.8356387 0.5492795
253239 0.9136306 -0.4065454
1871021 0.9973197 0.0731675
572511 -0.9832489 -0.1822682
831 0.1370879 0.9905589
301301 0.1638402 -0.9864869
29360 0.1433457 -0.9896727
28446 0.9008602 0.4341093
1679721 0.2425193 -0.9701466
649756 -0.8753416 -0.4835050
617123 0.8204468 -0.5717228
2109688 0.9379689 0.3467194
1981510 -0.9994219 0.0339979
2109687 0.2918990 0.9564491
35701 0.9029736 0.4296961
84032 0.0641110 -0.9979428
1510 0.0082848 -0.9999657
236753 -0.8810820 -0.4729636
1505 0.9134927 0.4068552
1496 -0.3762858 0.9265036
1511 0.9098324 0.4149760
2081703 -0.9827021 0.1851931
1736 0.3190652 -0.9477328
2184575 0.8277752 0.5610599
863 0.8249070 0.5652685
86170 0.8719464 -0.4896013
42838 0.5987104 0.8009656
1525 0.9199415 0.3920557
85874 0.8189921 0.5738048
1754 0.4147959 -0.9099145
129958 0.9173933 0.3979817
911092 0.4596580 0.8880960
499229 -0.0151625 -0.9998850
291990 0.9590783 0.2831409
1517 0.9866831 -0.1626542
44000 0.2319512 -0.9727274
252966 0.5417082 0.8405666
227388 0.8365948 -0.5478222
31909 -0.9806272 -0.1958835
28187 -0.8604307 -0.5095675
42422 -0.9512876 -0.3083050
375929 0.8600085 -0.5102797
69823 0.8541652 -0.5200017
158847 0.9897288 -0.1429577
484770 0.8810881 0.4729523
1702287 0.4263595 -0.9045538
2161821 -0.9974856 -0.0708696
1555112 0.9811835 0.1930775
1702221 0.9592456 -0.2825736
1712675 -0.9974934 0.0707597
1514105 0.4818243 0.8762678
39483 -0.2422601 0.9702113
1870984 0.9715718 -0.2367450
33033 -0.9995076 0.0313774
1260 -0.8445887 0.5354157
1912856 0.9279919 -0.3726003
1556 0.9338221 -0.3577377
1852374 -0.5061591 0.8624401
1871025 0.9905915 -0.1368519
1298 0.9988912 -0.0470776
332249 0.8907671 0.4544600
270 0.9091294 -0.4165138
1129 0.9548541 0.2970751
59930 0.9574278 0.2886729
1394889 0.9204486 0.3908636
1080068 0.9691042 -0.2466517
1209493 0.4998642 0.8661038
1173032 0.9705640 -0.2408434
82654 -0.0702293 0.9975309
1142 -0.2148534 0.9766463
2005458 0.8092474 0.5874680
142864 0.4147406 0.9099397
1186 0.9002746 0.4353224
373994 -0.3699040 0.9290700
1197 0.1369425 -0.9905790
1137095 0.9656807 0.2597323
136073 0.8975746 -0.4408625
109266 -0.9999567 0.0093101
312883 -0.9165076 0.4000174
1752063 0.8221857 0.5692193
482564 0.8979759 0.4400446
395961 0.9286268 0.3710152
1173027 0.9874130 -0.1581634
1206 0.9975443 -0.0700378
864702 0.8523935 -0.5229009
241425 -0.0277600 0.9996146
669359 0.9023334 0.4310388
2005460 0.4726560 0.8812470
1126 0.9375281 0.3479095
379064 0.9016513 0.4324639
65093 0.9897288 -0.1429577
713887 -0.3609746 -0.9325756
33072 0.9133245 0.4072325
54308 0.9369560 0.3494474
1807358 0.4928560 0.8701109
54299 0.8235486 0.5672458
1188228 0.9729792 0.2308928
2057 0.8352163 -0.5499216
500 0.3430239 0.9393267
120962 0.9754633 0.2201621
1108 0.8272343 0.5618571
1806508 0.9946992 -0.1028275
167964 0.9884409 -0.1516067
1986204 0.8948372 -0.4463926
133453 0.9927309 -0.1203551
1839801 0.9996413 -0.0267817
1005039 0.9584325 0.2853193
454171 -0.9991378 -0.0415175
2093 0.8771414 0.4802322
2130 0.3072179 -0.9516392
2132 0.8824197 0.4704630
2151 0.1361795 -0.9906842
214888 -0.9996168 0.0276822
38986 0.9859035 -0.1673148
1541959 -0.3810094 -0.9245712
35786 0.9908955 0.1346333
1636152 0.8006297 0.5991594
1891926 0.5762144 0.8172986
1331910 0.9650176 0.2621851
265606 -0.3031655 -0.9529379
120 0.9780453 -0.2083922
1387353 0.9294505 0.3689469
466153 0.8881154 0.4596206
128 -0.2503497 0.9681555
174633 0.9890211 0.1477743
1941349 0.8940103 0.4480464
107709 0.8563250 0.5164373
1838286 0.9029557 0.4297336
1796921 0.9458908 -0.3244851
395922 0.9876178 0.1568789
2026799 0.4839096 0.8751180
1704307 0.9908216 0.1351762
511746 0.9761035 -0.2173060
1307763 0.9908431 0.1350181
389348 -0.9715241 -0.2369409
83552 -0.9959124 0.0903241
71667 -0.8960445 0.4439642
83561 -0.8254627 -0.5644567
810 -0.9977357 0.0672563
1484118 0.9367910 0.3498893
1178516 0.5580741 0.8297912
651143 0.9989714 -0.0453447
94254 0.0721576 0.9973932
106 0.5501167 0.8350878
985 -0.9103525 -0.4138337
1784714 -0.3871397 0.9220211
316068 -0.9265817 -0.3760935
2183547 -0.8743213 -0.4853476
1795355 0.9879421 -0.1548239
1727163 0.9956053 -0.0936483
104 -0.4536991 -0.8911550
1191459 0.4941486 0.8693774
1006 -0.9811429 -0.1932838
1085624 0.2668203 0.9637463
247481 0.8156766 -0.5785082
281120 0.8512227 -0.5248047
29549 0.9901833 0.1397750
146919 0.9103357 0.4138708
2162713 -0.9694691 -0.2452135
253 -0.9536812 -0.3008193
1016 0.9124401 -0.4092103
252307 0.4200916 0.9074817
1486245 -0.2750302 -0.9614356
238 -0.9487111 0.3161444
754423 -0.1034435 0.9946353
1178778 -0.0711944 0.9974625
531844 -0.4645853 0.8855284
1249933 0.8879519 -0.4599362
1714848 -0.9814531 -0.1917022
256 -0.9514156 0.3079097
326319 -0.5940491 0.8044288
143223 0.4212952 -0.9069236
320324 -0.8000376 -0.5999498
983544 -0.9823412 0.1870982
2027857 -0.1605039 0.9870352
2282170 0.9948001 -0.1018464
101385 -0.5850729 0.8109807
1761453 0.9922806 0.1240128
1383885 0.9977360 -0.0672521
2058135 -0.9997061 0.0242436
2094025 -0.9136767 -0.4064417
616991 -0.1739186 -0.9847600
2069432 0.9211291 -0.3892573
63186 -0.8254627 -0.5644567
754429 -0.9530817 -0.3027132
1790137 -0.9897741 0.1426439
313588 -0.2868996 0.9579606
1453352 -0.8445887 0.5354157
1803846 -0.9930778 -0.1174586
1014 -0.9992663 0.0383008
762954 -0.1089770 0.9940443
1936081 0.1379897 -0.9904337
57029 -0.9861038 -0.1661306
191579 0.8910672 0.4538713
34098 0.2785773 -0.9604138
253245 -0.5219845 -0.8529550
336810 0.8019134 -0.5974403
1415657 0.8931981 0.4496633
242600 -0.9991704 0.0407241
1433126 -0.5777007 -0.8162487
52227 0.9584599 0.2852275
712710 -0.1542710 0.9880286
375288 0.8068637 0.5907376
1642646 -0.2231178 0.9747915
397865 0.5543871 0.8322590
1796646 -0.3470853 0.9378336
185300 0.9957231 -0.0923880
1400053 -0.8505540 0.5258877
1168034 0.4992043 0.8664843
1307839 -0.4650551 -0.8852817
354356 -0.9973627 0.0725789
2315862 -0.8961165 0.4438189
79329 -0.9365261 -0.3505980
446683 -0.9998444 0.0176386
1492898 -0.8866837 -0.4623764
1850526 -0.9998953 0.0144693
477680 -0.9999347 -0.0114240
2305508 -0.9999954 0.0030296
1986952 -0.3720285 -0.9282213
151895 0.9949347 -0.1005233
2350 0.8169096 -0.5767657
1457365 0.9816775 -0.1905498
274537 0.8617901 0.5072651
1096 0.9150673 -0.4033012
1100 0.9792749 0.2025357
281093 0.9493813 0.3141261
100716 0.9995887 -0.0286766
1134405 0.9994033 -0.0345392
861299 0.8484204 0.5293230
173480 0.9405960 0.3395279
658062 0.8895476 0.4568425
392734 -0.1437230 0.9896180
33075 0.8628124 0.5055243
458033 0.9018537 0.4320415
42253 0.9153643 0.4026266
180 0.9721127 0.2345141
28262 0.5088243 -0.8608704
154 0.8899182 0.4561202
81028 0.9773408 0.2116717
1131703 0.4335800 0.9011151
1307761 0.3544695 0.9350676
55206 0.9727462 -0.2318726
138 0.9863170 -0.1648600
64895 -0.1302277 -0.9914841
28452 0.9708945 0.2395077
29510 -0.2885923 -0.9574521
52584 0.8418189 -0.5397600
81412 0.9125850 -0.4088870
1197717 0.9200514 0.3917976
81468 0.4772261 -0.8787806
336261 0.0799315 -0.9968004
57487 0.2491098 -0.9684752
93466 0.9948891 -0.1009737
2420 0.9642507 0.2649918
1330330 0.8361818 0.5484523
1184387 0.5859955 0.8103143
160798 0.9999738 0.0072347
1006576 -0.2830481 0.9591057
69499 0.9971807 -0.0750377
381751 0.1342849 -0.9909428
63363 0.9475258 -0.3196792
168657 0.9955481 -0.0942551
940 -0.8445887 0.5354157
436114 0.9698017 -0.2438947
309805 -0.0711944 0.9974625
412593 0.9922640 -0.1241454
228745 0.9957749 -0.0918280
64160 0.5937962 0.8046155
851 -0.9368004 0.3498642
712368 -0.5050775 -0.8630740
826 0.9251074 -0.3797057
187101 -0.9009793 0.4338620
34105 -0.3810094 -0.9245712
936456 0.8646175 0.5024307
187145 0.8926146 0.4508205
423605 0.5655731 0.8246981
1408281 0.9129673 0.4080328
1295609 0.2941208 0.9557683
1653476 0.3670653 0.9301952
171695 0.9673074 -0.2536067
2352 0.8438626 -0.5365593
118000 0.4553540 -0.8903105
197162 -0.9168533 -0.3992244
513050 0.9485631 0.3165880
693075 0.8849511 0.4656840
634113 -0.9892963 -0.1459204
114186 -0.9987793 0.0493946
# Correr todo junto
# Definir parámetros
par(mar = c(5, 5, 1, 2) + 0.1)
# Inicializar el plot
plot(PCoA$points[ ,1], PCoA$points[ ,2], ylim = c(-0.5, 0.5),
xlab = paste("PCoA 1 (", explainedvar1, " %)", sep = ""),
ylab = paste("PCoA 2 (", explainedvar2, " %)", sep = ""),
pch = 5, cex = 1.0, type = "n", cex.lab = 1.0, cex.axis = 1.2, axes = FALSE)
# Añadir ejes
axis(side = 1, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
axis(side = 2, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
abline(h = 0, v = 0, lty = 3)
box(lwd = 2)
# Añadir puntos y etiquetas
points(PCoA$points[ ,1], PCoA$points[ ,2],
pch = 19, cex = 3, bg = "blue", col = "blue")
text(PCoA$points[ ,1], PCoA$points[ ,2],
labels = row.names(PCoA$points))
CalakmulREL <- abund_table
for(i in 1:nrow(abund_table)){
CalakmulREL[i, ] = abund_table[i, ] / sum(abund_table[i, ])
}
# Calcular y añadir el score de las especies
PCoA <- add.spec.scores(PCoA,CalakmulREL,method = "pcoa.scores",Rscale=TRUE,scaling=1, multi=1)
text(PCoA$cproj[ ,1], PCoA $cproj[ ,2],
labels = row.names(PCoA$cproj),cex=0.5, col = "blue")

genus_corr <- add.spec.scores(PCoA, CalakmulREL, method = "cor.scores")$cproj
corrcut <- 0.7 # definir el cutoff
import_genus <- genus_corr[abs(genus_corr[, 1]) >= corrcut | abs(genus_corr[, 2]) >= corrcut, ]
# Código para los generos con un p-valor<0.05
fit <- envfit(PCoA, CalakmulREL, perm = 999)
## 'nperm' >= set of all permutations: complete enumeration.
## Set of permutations < 'minperm'. Generating entire set.
## 'nperm' >= set of all permutations: complete enumeration.
## Set of permutations < 'minperm'. Generating entire set.
plot(fit, p.max = 0.05, cex=0.5, col = "red")

Top10

bc_dist <- vegdist(Abund10, "bray")
PCoA <- cmdscale(bc_dist, eig = TRUE, k = 2)
PCoA
## $points
##               [,1]         [,2]
## 374     0.44387444  0.039250552
## 407    -0.20864520 -0.105913581
## 32008  -0.14811785 -0.006984139
## 286    -0.03651816  0.098916814
## 1883    0.55050030 -0.091453429
## 1763   -0.08427325  0.062764785
## 36814  -0.12542113  0.021382746
## 1873   -0.02717031  0.104829143
## 196162 -0.20621583 -0.099970376
## 191495 -0.15801301 -0.022822515
## 
## $eig
##  [1]  6.579441e-01  5.685607e-02  4.621780e-03  4.552168e-03  3.788089e-03
##  [6]  5.087640e-04  2.481302e-04 -6.938894e-18 -5.051696e-05 -3.836861e-03
## 
## $x
## NULL
## 
## $ac
## [1] 0
## 
## $GOF
## [1] 0.9759610 0.9811687
# Variación explicada primer eje
explainedvar1 <- round(PCoA$eig[1] / sum(PCoA$eig), 2) * 100
explainedvar1
## [1] 91
# Variación explicada segundo eje
explainedvar2 <- round(PCoA$eig[2] / sum(PCoA$eig), 2) * 100
explainedvar2
## [1] 8
# Suma
sum_eig <- sum(explainedvar1, explainedvar2)
sum_eig
## [1] 99

Criterios de evaluación

# Correr todo junto
# Definir los parametros del plot
par(mar = c(5, 5, 1, 2) + 0.1)
# Fraficar eigenvalores
plot(PCoA$eig, xlab = "PCoA", ylab = "Eigenvalue",
las = 1, cex.lab = 1.5, pch = 16)
# Añadir expextación del criterio de Kaiser-Guttman y el modelo Broken Stick
abline(h = mean(PCoA$eig), lty = 2, lwd = 2, col = "blue")
b_stick <- bstick(8, sum(PCoA$eig))
lines(1:8, b_stick, type = "l", lty = 4, lwd = 2, col = "red")
# Añadir legendas
legend("topright", legend = c("Avg Eigenvalue", "Broken-Stick"),
lty = c(2, 4), bty = "n", col = c("blue", "red"))

Gráfico

# Correr todo junto
# Definir parámetros
par(mar = c(5, 5, 1, 2) + 0.1)
# Initiate Plot
plot(PCoA$points[ ,1], PCoA$points[ ,2], ylim = c(-0.5, 0.5),
xlab = paste("PCoA 1 (", explainedvar1, " %)", sep = ""),
ylab = paste("PCoA 2 (", explainedvar2, " %)", sep = ""),
pch = 5, cex = 1.0, type = "n", cex.lab = 1.0, cex.axis = 1.2, axes = FALSE)
# Añadir ejes
axis(side = 1, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
axis(side = 2, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
abline(h = 0, v = 0, lty = 3)
box(lwd = 2)
# Añadir puntos y texto
points(PCoA$points[ ,1], PCoA$points[ ,2],
pch = 19, cex = 3, bg = "blue", col = "blue")
text(PCoA$points[ ,1], PCoA$points[ ,2],
labels = row.names(PCoA$points))

# Correr todo junto
# Definir parámetros
par(mar = c(5, 5, 1, 2) + 0.1)
# Inicializar el plot
plot(PCoA$points[ ,1], PCoA$points[ ,2], ylim = c(-0.5, 0.5),
xlab = paste("PCoA 1 (", explainedvar1, " %)", sep = ""),
ylab = paste("PCoA 2 (", explainedvar2, " %)", sep = ""),
pch = 5, cex = 1.0, type = "n", cex.lab = 1.0, cex.axis = 1.2, axes = FALSE)
# Añadir ejes
axis(side = 1, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
axis(side = 2, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
abline(h = 0, v = 0, lty = 3)
box(lwd = 2)
# Añadir puntos y etiquetas
points(PCoA$points[ ,1], PCoA$points[ ,2],
pch = 19, cex = 3, bg = "blue", col = "blue")
text(PCoA$points[ ,1], PCoA$points[ ,2],
labels = row.names(PCoA$points))
CalakmulREL <- Abund10
for(i in 1:nrow(Abund10)){
CalakmulREL[i, ] = Abund10[i, ] / sum(Abund10[i, ])
}
# Calcular y añadir el score de las especies
PCoA <- add.spec.scores(PCoA,CalakmulREL,method = "pcoa.scores",Rscale=TRUE,scaling=1, multi=1)
text(PCoA$cproj[ ,1], PCoA $cproj[ ,2],
labels = row.names(PCoA$cproj),cex=0.5, col = "blue")

genus_corr <- add.spec.scores(PCoA, CalakmulREL, method = "cor.scores")$cproj
corrcut <- 0.8 #definir el cutoff  
import_genus <- genus_corr[abs(genus_corr[, 1]) >= corrcut | abs(genus_corr[, 2]) >= corrcut, ]
# Los 12 géneros importantes con correlación mayor o igual a 0,8 
# a lo largo de los ejes PCoA se imprimen a continuación:
import_genus[complete.cases(import_genus),] %>%
kbl(booktabs = TRUE, align = "c",
caption = "Géneros más importantes con correlación mayor o igual a 0.8" ) %>%
kable_styling(position = "center",
latex_options = c("hold_position", "striped"),
font_size = 7)
Géneros más importantes con correlación mayor o igual a 0.8
Dim1 Dim2
# Correr todo junto
# Definir parámetros
par(mar = c(5, 5, 1, 2) + 0.1)
# Inicializar el plot
plot(PCoA$points[ ,1], PCoA$points[ ,2], ylim = c(-0.5, 0.5),
xlab = paste("PCoA 1 (", explainedvar1, " %)", sep = ""),
ylab = paste("PCoA 2 (", explainedvar2, " %)", sep = ""),
pch = 5, cex = 1.0, type = "n", cex.lab = 1.0, cex.axis = 1.2, axes = FALSE)
# Añadir ejes
axis(side = 1, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
axis(side = 2, labels = T, lwd.ticks = 2, cex.axis = 1.2, las = 1)
abline(h = 0, v = 0, lty = 3)
box(lwd = 2)
# Añadir puntos y etiquetas
points(PCoA$points[ ,1], PCoA$points[ ,2],
pch = 19, cex = 2, bg = "blue", col = "blue")
text(PCoA$points[ ,1], PCoA$points[ ,2],
labels = row.names(PCoA$points))
CalakmulREL <- Abund10
for(i in 1:nrow(Abund10)){
CalakmulREL[i, ] = Abund10[i, ] / sum(Abund10[i, ])
}
# Calcular y añadir el score de las especies
PCoA <- add.spec.scores(PCoA,CalakmulREL,method = "pcoa.scores",Rscale=TRUE,scaling=1, multi=1)
text(PCoA$cproj[ ,1], PCoA $cproj[ ,2],
labels = row.names(PCoA$cproj),cex=0.5, col = "blue")

genus_corr <- add.spec.scores(PCoA, CalakmulREL, method = "cor.scores")$cproj
corrcut <- 0.7 # definir el cutoff
import_genus <- genus_corr[abs(genus_corr[, 1]) >= corrcut | abs(genus_corr[, 2]) >= corrcut, ]
# Código para los generos con un p-valor<0.05
fit <- envfit(PCoA, CalakmulREL, perm = 999)
plot(fit, p.max = 0.05, cex=0.5, col = "red")

Escalamiento multidimensional no métrico (NMDS)

Todos

bc_nmds <- metaMDS(abund_table, dist = "bray")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0 
## Run 1 stress 0 
## ... Procrustes: rmse 0.06591674  max resid 0.07489395 
## Run 2 stress 0 
## ... Procrustes: rmse 0.2926706  max resid 0.4058441 
## Run 3 stress 0 
## ... Procrustes: rmse 0.03262501  max resid 0.03866392 
## Run 4 stress 0 
## ... Procrustes: rmse 0.01431391  max resid 0.0161337 
## Run 5 stress 0 
## ... Procrustes: rmse 0.05551793  max resid 0.05922182 
## Run 6 stress 0 
## ... Procrustes: rmse 0.2108313  max resid 0.2880816 
## Run 7 stress 0 
## ... Procrustes: rmse 0.2702235  max resid 0.3552578 
## Run 8 stress 0 
## ... Procrustes: rmse 0.3086159  max resid 0.4006153 
## Run 9 stress 0 
## ... Procrustes: rmse 0.033168  max resid 0.03667409 
## Run 10 stress 0 
## ... Procrustes: rmse 0.295754  max resid 0.4165082 
## Run 11 stress 0 
## ... Procrustes: rmse 0.03418084  max resid 0.03841165 
## Run 12 stress 0 
## ... Procrustes: rmse 0.07220893  max resid 0.0791087 
## Run 13 stress 0 
## ... Procrustes: rmse 0.2715364  max resid 0.3538363 
## Run 14 stress 0 
## ... Procrustes: rmse 0.09322069  max resid 0.1158608 
## Run 15 stress 0 
## ... Procrustes: rmse 0.230417  max resid 0.3173905 
## Run 16 stress 0 
## ... Procrustes: rmse 0.3036194  max resid 0.3854338 
## Run 17 stress 0 
## ... Procrustes: rmse 0.06612041  max resid 0.07854557 
## Run 18 stress 0 
## ... Procrustes: rmse 0.2759622  max resid 0.3811386 
## Run 19 stress 0 
## ... Procrustes: rmse 0.1126835  max resid 0.1425449 
## Run 20 stress 0 
## ... Procrustes: rmse 0.1246335  max resid 0.1531633 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress < smin
## Warning in metaMDS(abund_table, dist = "bray"): stress is (nearly) zero: you may
## have insufficient data
## Warning in postMDS(out$points, dis, plot = max(0, plot - 1), ...): skipping
## half-change scaling: too few points below threshold
bc_nmds
## 
## Call:
## metaMDS(comm = abund_table, distance = "bray") 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(abund_table)) 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0 
## Stress type 1, weak ties
## Best solution was not repeated after 20 tries
## The best solution was from try 0 (metric scaling or null solution)
## Scaling: centring, PC rotation 
## Species: expanded scores based on 'wisconsin(sqrt(abund_table))'
# Graficamos
ordiplot (bc_nmds, type = 't')

ordiplot(bc_nmds, display = "sites", type = "text")

par (mfrow = c(1,2))
stressplot (bc_nmds)
plot (bc_nmds, display = 'sites', type = 't', main = 'Goodness of fit')
points (bc_nmds, display = 'sites', cex = goodness (bc_nmds)*300)

Top10

bc_nmds <- metaMDS(Abund10, dist = "bray")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0 
## Run 1 stress 0.00009858673 
## ... Procrustes: rmse 0.06846217  max resid 0.09471042 
## Run 2 stress 0.00007644814 
## ... Procrustes: rmse 0.02279074  max resid 0.03422784 
## Run 3 stress 0.00009120155 
## ... Procrustes: rmse 0.06435496  max resid 0.08031878 
## Run 4 stress 0.00007350407 
## ... Procrustes: rmse 0.04773922  max resid 0.07144841 
## Run 5 stress 0.00009136378 
## ... Procrustes: rmse 0.02906196  max resid 0.04328228 
## Run 6 stress 0.00009974308 
## ... Procrustes: rmse 0.06806319  max resid 0.09774495 
## Run 7 stress 0.04472006 
## Run 8 stress 0.00009132276 
## ... Procrustes: rmse 0.05424203  max resid 0.08680448 
## Run 9 stress 0.0000971888 
## ... Procrustes: rmse 0.06462366  max resid 0.1077245 
## Run 10 stress 0.2531766 
## Run 11 stress 0.07322694 
## Run 12 stress 0.00009842633 
## ... Procrustes: rmse 0.06636149  max resid 0.1127756 
## Run 13 stress 0.00008811215 
## ... Procrustes: rmse 0.06832711  max resid 0.08582095 
## Run 14 stress 0.00006683377 
## ... Procrustes: rmse 0.0612015  max resid 0.09595884 
## Run 15 stress 0.00007242056 
## ... Procrustes: rmse 0.04876336  max resid 0.06714538 
## Run 16 stress 0.00009006661 
## ... Procrustes: rmse 0.02525963  max resid 0.04679179 
## Run 17 stress 0.00009070948 
## ... Procrustes: rmse 0.04628716  max resid 0.07633794 
## Run 18 stress 0.00009925197 
## ... Procrustes: rmse 0.0682848  max resid 0.09791053 
## Run 19 stress 0.0000962134 
## ... Procrustes: rmse 0.02523965  max resid 0.03953934 
## Run 20 stress 0.00008270894 
## ... Procrustes: rmse 0.02414048  max resid 0.03489533 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     17: stress < smin
##      2: stress ratio > sratmax
##      1: scale factor of the gradient < sfgrmin
## Warning in metaMDS(Abund10, dist = "bray"): stress is (nearly) zero: you may
## have insufficient data
bc_nmds
## 
## Call:
## metaMDS(comm = Abund10, distance = "bray") 
## 
## global Multidimensional Scaling using monoMDS
## 
## Data:     wisconsin(sqrt(Abund10)) 
## Distance: bray 
## 
## Dimensions: 2 
## Stress:     0 
## Stress type 1, weak ties
## Best solution was not repeated after 20 tries
## The best solution was from try 0 (metric scaling or null solution)
## Scaling: centring, PC rotation, halfchange scaling 
## Species: expanded scores based on 'wisconsin(sqrt(Abund10))'
# Graficamos
ordiplot (bc_nmds, type = 't')

ordiplot(bc_nmds, display = "sites", type = "text")

par (mfrow = c(1,2))
stressplot (bc_nmds)
plot (bc_nmds, display = 'sites', type = 't', main = 'Goodness of fit')
points (bc_nmds, display = 'sites', cex = goodness (bc_nmds)*300)

Análisis de correspondencia

Todos

Calakmul_genus_cca <- cca(abund_table)
Calakmul_genus_cca
## Call: cca(X = abund_table)
## 
##               Inertia Rank
## Total         0.00947     
## Unconstrained 0.00947    2
## Inertia is scaled Chi-square 
## 
## Eigenvalues for unconstrained axes:
##      CA1      CA2 
## 0.008616 0.000854
plot(Calakmul_genus_cca, display="sites")

plot(Calakmul_genus_cca, display="sites", type="p")

ordiplot(Calakmul_genus_cca)

evplot <- function(ev)
{# Broken stick model (MacArthur 1957)
n <- length(ev)
bsm <- data.frame(j=seq(1:n), p=0)
bsm$p[1] <- 1/n
for (i in 2:n) bsm$p[i] <- bsm$p[i-1] + (1/(n + 1 - i))
bsm$p <- 100*bsm$p/n
# Graficar eigenvalores y porcentaje de varianza en cada eje
op <- par(mfrow=c(2,1))
barplot(ev, main="Eigenvalues", col="bisque", las=2)
abline(h=mean(ev), col="red")
legend("topright", "Average eigenvalue", lwd=1, col=2, bty="n")
barplot(t(cbind(100*ev/sum(ev), bsm$p[n:1])), beside=TRUE,
main=" % variation", col=c("bisque",2), las=2)
legend("topright", c(" % eigenvalue", "Broken stick model"),
pch=15, col=c("bisque",2), bty="n")
par(op)
}
# Plot de eigenvalores y porcentaje de varianza en cada eje
ev <- Calakmul_genus_cca$CA$eig

evplot(ev)

Top10

Calakmul_genus_cca <- cca(Abund10)
Calakmul_genus_cca
## Call: cca(X = Abund10)
## 
##                Inertia Rank
## Total         0.003977     
## Unconstrained 0.003977    2
## Inertia is scaled Chi-square 
## 
## Eigenvalues for unconstrained axes:
##      CA1      CA2 
## 0.003564 0.000414
plot(Calakmul_genus_cca, display="sites")

plot(Calakmul_genus_cca, display="sites", type="p")

ordiplot(Calakmul_genus_cca, type = "text") 

ordiplot(Calakmul_genus_cca) 

evplot <- function(ev)
{# Broken stick model (MacArthur 1957)
n <- length(ev)
bsm <- data.frame(j=seq(1:n), p=0)
bsm$p[1] <- 1/n
for (i in 2:n) bsm$p[i] <- bsm$p[i-1] + (1/(n + 1 - i))
bsm$p <- 100*bsm$p/n
# Graficar eigenvalores y porcentaje de varianza en cada eje
op <- par(mfrow=c(2,1))
barplot(ev, main="Eigenvalues", col="bisque", las=2)
abline(h=mean(ev), col="red")
legend("topright", "Average eigenvalue", lwd=1, col=2, bty="n")
barplot(t(cbind(100*ev/sum(ev), bsm$p[n:1])), beside=TRUE,
main=" % variation", col=c("bisque",2), las=2)
legend("topright", c(" % eigenvalue", "Broken stick model"),
pch=15, col=c("bisque",2), bty="n")
par(op)
}
# Plot de eigenvalores y porcentaje de varianza en cada eje
ev <- Calakmul_genus_cca$CA$eig

evplot(ev)